A method for automatically detecting control model parameters of a closed-loop numerical control system in laser cutting

Abstract

This invention relates to the field of laser cutting control technology, and more particularly to an automatic detection method for control model parameters of a closed-loop CNC system in laser cutting. The method includes the following steps: acquiring a raw system state dataset; extracting interference features based on the raw system state dataset to obtain system interference feature data, wherein the system interference feature data includes steam plume interference feature data and servo motor thermal feature data; extracting interference features from the steam plume interference feature data to obtain plume interference frequency domain characteristic data; performing motor parameter drift analysis on the servo motor thermal characteristic data to obtain motor parameter drift characteristic data; and constructing a system parameter drift characteristic library based on the plume interference frequency domain characteristic data and the motor parameter drift characteristic data. This invention effectively solves the cutting accuracy problem caused by parameter drift and significantly improves the overall performance and reliability of the closed-loop CNC system in laser cutting.

Description

Technical Field

[0001] This invention relates to the field of laser cutting control technology, and in particular to an automatic detection method for control model parameters of a closed-loop CNC system in laser cutting. Background Technology

[0002] In modern manufacturing, laser cutting technology, as one of the advanced processing methods, is widely used in the manufacturing of automotive exhaust systems, especially in the high-precision contour machining of stainless steel exhaust manifolds. With the continuous development of the automotive industry, the requirements for the machining precision of exhaust system components are increasing, and laser cutting technology, with its high precision, high efficiency, and good adaptability, has become an ideal choice for machining stainless steel exhaust manifolds.

[0003] Stainless steel exhaust manifolds typically feature complex spatial curved flange interfaces and thin-walled corrugated pipe sections. Their assembly airtightness requires the cut kerf width error to be stably controlled within ±0.03mm. However, during long-term continuous machining of high-temperature resistant stainless steel (such as 409L) with titanium alloy coatings, the dynamic response parameters of the closed-loop CNC system control model in laser cutting exhibit progressive drift. This drift primarily stems from two rarely monitored coupling factors: first, the intermittent optical interference of the metal vapor plume generated by the high-temperature laser ablation of the coating on the optical encoder reading head, resulting in millisecond-level micro-jitter in the position feedback signal; second, the attenuation of the electromagnetic constant of the servo motor due to the accumulation of winding temperature during frequent starts and stops of small-curvature arcs. Summary of the Invention

[0004] Based on this, it is necessary for the present invention to provide an automatic detection method for control model parameters of a closed-loop CNC system in laser cutting, so as to solve at least one of the above-mentioned technical problems.

[0005] To achieve the above objectives, an automatic detection method for control model parameters of a closed-loop CNC system in laser cutting is provided, comprising the following steps:

[0006] Step S1: Obtain the original system state dataset; extract disturbance features based on the original system state dataset to obtain system disturbance feature data, which includes steam plume disturbance feature data and servo motor thermal feature data;

[0007] Step S2: Extract interference features from the steam plume interference characteristic data to obtain plume interference frequency domain characteristic data; perform motor parameter drift analysis on the servo motor thermal characteristic data to obtain motor parameter drift characteristic data; construct a system parameter drift characteristic library based on the plume interference frequency domain characteristic data and the motor parameter drift characteristic data.

[0008] Step S3: Based on the system parameter drift symptom library, perform online identification and weight adjustment of the system feedforward compensation coefficients to obtain weighted CNC system feedforward compensation data; perform parameter oscillation suppression on the weighted CNC system feedforward compensation data to obtain stable feedforward compensation parameters;

[0009] Step S4: Quantify the system stability margin based on the stable feedforward compensation parameters to obtain the stability margin of the closed-loop CNC system; perform iterative optimization of the speed loop gain parameters based on the stability margin of the closed-loop CNC system to obtain the optimized speed loop control parameters.

[0010] Step S5: Predict the cutting speed based on the speed loop optimization control parameters to obtain the predicted cutting speed distribution data; perform predictive compensation control on the current cutting trajectory data based on the predicted cutting speed distribution data to obtain the final cutting trajectory control data.

[0011] This invention solves the traditional monitoring challenges by collecting and analyzing raw system state datasets and acquiring high-precision raw data using multiple sensors. It accurately identifies steam plume interference characteristics and servo motor thermal characteristics. By online identification and adjustment of the system's feedforward compensation coefficients combined with parameter oscillation suppression, the system's adaptability and stability to parameter drift are effectively improved, ensuring a smooth cutting process. Quantifying the system's stability margin and iteratively optimizing the speed loop gain parameters enhances the system's dynamic response performance, achieving efficient and stable cutting control. Predictive compensation control based on cutting speed and current cutting trajectory data significantly improves the accuracy and precision of the cutting trajectory, meeting the requirements of high-precision cutting. This invention effectively solves the cutting accuracy problem caused by parameter drift, significantly improves the overall performance and reliability of closed-loop CNC systems in laser cutting, reduces the defect rate, and improves production efficiency and product quality.

[0012] Preferably, step S1 includes the following steps:

[0013] Step S11: Initialize the sensor network configuration of the closed-loop CNC system in laser cutting to obtain the sensor network configuration parameters. The sensor network initialization configuration includes setting the sampling frequency of the optical encoder to 10kHz, the sampling frequency of the servo motor current sensor to 5kHz, and the sampling frequency of the vibration accelerometer to 2kHz. At the same time, configure the CAN bus communication protocol baud rate to 1Mbps.

[0014] Step S12: Acquire position signals from the optical encoder based on the sensor network configuration parameters to obtain the raw position feedback signal;

[0015] Step S13: Monitor the three-phase current of the servo motor based on the original position feedback signal to obtain the characteristic data of the motor drive current;

[0016] Step S14: Based on the characteristic data of the motor drive current, perform triaxial monitoring of the vibration acceleration of the cutting head to obtain the dynamic response data of the cutting head;

[0017] Step S15: Based on the dynamic response data of the cutting head, monitor the environmental parameters of the closed-loop CNC system in laser cutting in real time to obtain environmental impact factor data;

[0018] Step S16: Based on the environmental impact factor data, timestamp the signals of each sensor to obtain time-stamped multi-source sensor data, where each sensor signal includes the optical encoder position signal, the servo motor three-phase current signal and the cutting head vibration acceleration signal.

[0019] Step S17: Generate the original system state dataset based on time-stamped multi-source sensor data;

[0020] Step S18: Extract disturbance features based on the original system state dataset to obtain system disturbance feature data, which includes steam plume disturbance feature data and servo motor thermal feature data.

[0021] Preferably, step S18 includes the following steps:

[0022] Step S181: Perform time series interpolation on the original system state dataset to obtain a system state dataset with a uniform sampling rate;

[0023] Step S182: Perform time window segmentation based on the unified sampling rate system state dataset to obtain synchronized system state data;

[0024] Step S183: Initialize the fiber optic sensor with wavelength demodulation based on the synchronization system status data to obtain the reference wavelength of the fiber optic sensor; perform spectral feature detection on the metal vapor plume based on the reference wavelength of the fiber optic sensor to obtain the plume spectral absorption data;

[0025] Step S184: Perform vapor density inversion based on plume spectral absorption data to obtain vapor plume interference characteristic data;

[0026] Step S185: Reconstruct the temperature field of the thermocouple array based on the synchronous system state data to obtain the servo motor temperature field distribution data;

[0027] Step S186: Calculate the temperature gradient of the motor windings based on the servo motor temperature field distribution data to obtain the servo motor thermal characteristic data.

[0028] Preferably, step S2, which involves extracting interference features from the steam plume interference data, includes:

[0029] Frequency domain transformation preprocessing is performed based on steam plume interference characteristic data to obtain plume frequency domain transformation preprocessing data. The frequency domain transformation preprocessing includes windowing the steam plume interference characteristic data using a Hamming window function, with the window length set to 1024 sampling points and an overlap rate of 75%. The data length is extended to 2048 points using zero-filling technology.

[0030] Fast Fourier Transform is performed on the preprocessed data of the plume frequency domain transform to obtain the frequency domain data of the position feedback signal;

[0031] Amplitude and phase spectra are separated based on the frequency domain data of the position feedback signal to obtain separated frequency domain feature data;

[0032] Power spectral density is estimated based on the separated frequency domain characteristic data to obtain the power spectral data of the position feedback signal;

[0033] Based on the power spectrum data of the location feedback signal, spectral peak detection and location are performed to obtain frequency domain characteristic data of plume interference.

[0034] Preferably, step S2, which involves analyzing the motor parameter drift of the servo motor thermal characteristic data, includes:

[0035] Wavelet basis function selection is performed based on servo motor thermal characteristic data to obtain wavelet transform parameters. The wavelet basis function selection includes selecting Daubechies4 wavelet as mother wavelet function, setting the number of decomposition layers to 6, the scale parameter a varying from 1 to 64, and the step size of the translation parameter b set to 1 / 10 of the sampling interval.

[0036] Continuous wavelet transform is performed on the thermal characteristic data of the servo motor according to the wavelet transform parameters to obtain the time-frequency domain data of the motor temperature.

[0037] Electromagnetic constant attenuation is identified based on motor temperature time-frequency domain data to obtain raw data of electromagnetic constant attenuation.

[0038] Motor thermal characteristics are identified based on the original data of electromagnetic constant decay to obtain instantaneous parameters of motor thermal characteristics; Huang transform decomposition is performed on the instantaneous parameters of motor thermal characteristics to obtain the intrinsic mode function data of thermal characteristics.

[0039] Instantaneous frequency calculations are performed based on thermal characteristic intrinsic mode function data to obtain motor parameter drift characteristic data.

[0040] Preferably, the construction of the system parameter drift indicator library based on the plume interference frequency domain indicator data and the motor parameter drift indicator data in step S2 includes:

[0041] Construct a plume interference feature vector based on plume interference frequency domain characteristic data;

[0042] Construct a motor drift feature vector based on motor parameter drift symptom data;

[0043] Multidimensional feature space mapping is performed based on the plume interference feature vector and the motor drift feature vector to obtain the fused system state feature space data;

[0044] Clustering and pattern recognition are performed on the fusion system state feature space data to obtain a system parameter drift sign library.

[0045] Preferably, step S3 includes the following steps:

[0046] Step S31: Initialize the system parameters using the recursive least squares method based on the system parameter drift symptom library to obtain the initial parameters of the recursive algorithm;

[0047] Step S32: Real-time data acquisition is performed on the closed-loop CNC system in laser cutting to obtain system input and output data; regression matrix is ​​constructed on the system input and output data based on the initial parameters of the recursive algorithm to obtain the regression matrix of the closed-loop CNC system;

[0048] Step S33: Perform recursive least squares parameter updates based on the regression matrix of the closed-loop CNC system to obtain real-time CNC system parameter estimation data;

[0049] Step S34: Extract the feedforward compensation coefficients based on the parameter estimation data of the real-time CNC system to obtain the original data of the feedforward compensation coefficients;

[0050] Step S35: Calculate the forgetting factor based on the original data of the feedforward compensation coefficient to obtain the forgetting factor data of the CNC system;

[0051] Step S36: Adjust the weights of the original data of the feedforward compensation coefficients based on the forgetting factor data of the CNC system to obtain the weighted feedforward compensation data of the CNC system;

[0052] Step S37: Perform parameter oscillation suppression on the feedforward compensation data of the weighted CNC system to obtain stable feedforward compensation parameters.

[0053] Preferably, step S37 includes the following steps:

[0054] Step S371: Extract historical data from the sliding window based on the weighted CNC system feedforward compensation data to obtain a set of historical CNC system data windows;

[0055] Step S372: Perform weighted least squares fitting based on the historical CNC system data window set to obtain the feedforward compensation coefficients of the CNC system;

[0056] Step S373: Calculate the parameter change rate based on the feedforward compensation coefficient of the CNC system to obtain the compensation coefficient change rate data;

[0057] Step S374: Perform threshold detection and amplitude limiting based on the rate of change of the compensation coefficient to obtain the corrected feedforward compensation coefficient;

[0058] Step S375: Perform system stability assessment based on the modified feedforward compensation coefficient to obtain system stability assessment data; select an adjustment strategy based on the system stability assessment data and the preset adjustment strategy set to obtain the effective strategy adjustment parameters of the system;

[0059] Step S376: Adjust the parameters based on the effective strategy of the system to perform parameter convergence verification and obtain the convergence verification data of the closed-loop CNC system;

[0060] Step S377: Confirm the final parameters based on the convergence verification data of the closed-loop CNC system to obtain stable feedforward compensation parameters.

[0061] Preferably, step S4 includes the following steps:

[0062] Step S41: Construct a mathematical model of the closed-loop system based on the stable feedforward compensation parameters to obtain the system state-space model;

[0063] Step S42: Extract the characteristic polynomial coefficients based on the system state-space model to obtain the characteristic polynomial parameters of the closed-loop system;

[0064] Step S43: Identify the system pole locations based on the characteristic polynomial parameters of the closed-loop system to obtain the closed-loop system pole data;

[0065] Step S44: Identify stability criteria based on the closed-loop system pole data to obtain the stability coefficients of the closed-loop system;

[0066] Step S45: Evaluate the stability margin of the closed-loop CNC system in laser cutting based on the stability coefficient of the closed-loop system to obtain the stability margin of the closed-loop CNC system.

[0067] Step S46: Iteratively optimize the speed loop gain parameters based on the stability margin of the closed-loop CNC system to obtain the optimized speed loop control parameters.

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

[0069] Step S51: Acquire the original trajectory data of the closed-loop CNC system in laser cutting to obtain the original cutting trajectory data; Based on the speed loop optimization control parameters, reconstruct the original cutting trajectory data to obtain the standard cutting trajectory data.

[0070] Step S52: Calculate the curvature based on the standard cutting trajectory data to obtain the cutting trajectory curvature distribution data;

[0071] Step S53: Extract geometric feature parameters based on the curvature distribution data of the cutting trajectory to obtain the geometric feature data of the cutting trajectory;

[0072] Step S54: Calculate the rate of change of tangent direction based on the geometric feature data of the cutting trajectory to obtain the feature data of the change of cutting direction;

[0073] Step S55: Perform velocity planning and prediction based on the cutting direction change feature data to obtain the cutting predicted velocity distribution data;

[0074] Step S56: Real-time trajectory data acquisition is performed on the closed-loop CNC system during laser cutting to obtain the current cutting trajectory data;

[0075] Step S57: Perform predictive compensation control on the current cutting trajectory data based on the cutting prediction speed distribution data to obtain the final cutting trajectory control data. Attached Figure Description

[0076] Other features, objects, and advantages of the present invention will become more apparent from the following detailed description taken in conjunction with the accompanying drawings:

[0077] Figure 1 A flowchart illustrating the steps of an automatic detection method for control model parameters of a closed-loop CNC system in laser cutting is shown in one embodiment.

[0078] Figure 2 A detailed flowchart of step S3 of one embodiment is shown;

[0079] Figure 3 A detailed flowchart of step S37 of one embodiment is shown. Detailed Implementation

[0080] The technical method of the present invention will now be clearly and completely described with reference to the accompanying drawings. Obviously, the described embodiments are only some, not all, of the embodiments of the present invention. All other embodiments obtained by those skilled in the art based on the embodiments of the present invention without inventive effort are within the scope of protection of the present invention.

[0081] Furthermore, the accompanying drawings are merely illustrative of the invention and are not necessarily drawn to scale. The same reference numerals in the drawings denote the same or similar parts, and therefore repeated descriptions of them will be omitted. Some block diagrams shown in the drawings are functional entities and do not necessarily correspond to physically or logically independent entities. These functional entities can be implemented in software, in one or more hardware modules or integrated circuits, or in different network and / or processor methods and / or microcontroller methods.

[0082] It should be understood that although the terms "first," "second," etc., may be used herein to describe various units, these units should not be limited by these terms. These terms are used merely to distinguish one unit from another. For example, without departing from the scope of the exemplary embodiments, a first unit may be referred to as a second unit, and similarly, a second unit may be referred to as a first unit. The term "and / or" as used herein includes any and all combinations of one or more of the associated listed items.

[0083] To achieve the above objectives, please refer to Figures 1 to 3 This invention provides an automatic detection method for control model parameters of a closed-loop CNC system in laser cutting, comprising the following steps:

[0084] Step S1: Obtain the original system state dataset; extract disturbance features based on the original system state dataset to obtain system disturbance feature data, which includes steam plume disturbance feature data and servo motor thermal feature data;

[0085] Step S2: Extract interference features from the steam plume interference characteristic data to obtain plume interference frequency domain characteristic data; perform motor parameter drift analysis on the servo motor thermal characteristic data to obtain motor parameter drift characteristic data; construct a system parameter drift characteristic library based on the plume interference frequency domain characteristic data and the motor parameter drift characteristic data.

[0086] Step S3: Based on the system parameter drift symptom library, perform online identification and weight adjustment of the system feedforward compensation coefficients to obtain weighted CNC system feedforward compensation data; perform parameter oscillation suppression on the weighted CNC system feedforward compensation data to obtain stable feedforward compensation parameters;

[0087] Step S4: Quantify the system stability margin based on the stable feedforward compensation parameters to obtain the stability margin of the closed-loop CNC system; perform iterative optimization of the speed loop gain parameters based on the stability margin of the closed-loop CNC system to obtain the optimized speed loop control parameters.

[0088] Step S5: Predict the cutting speed based on the speed loop optimization control parameters to obtain the predicted cutting speed distribution data; perform predictive compensation control on the current cutting trajectory data based on the predicted cutting speed distribution data to obtain the final cutting trajectory control data.

[0089] Preferably, step S1 includes the following steps:

[0090] Step S11: Initialize the sensor network configuration of the closed-loop CNC system in laser cutting to obtain the sensor network configuration parameters. The sensor network initialization configuration includes setting the sampling frequency of the optical encoder to 10kHz, the sampling frequency of the servo motor current sensor to 5kHz, and the sampling frequency of the vibration accelerometer to 2kHz. At the same time, configure the CAN bus communication protocol baud rate to 1Mbps.

[0091] Step S12: Acquire position signals from the optical encoder based on the sensor network configuration parameters to obtain the raw position feedback signal;

[0092] Step S13: Monitor the three-phase current of the servo motor based on the original position feedback signal to obtain the characteristic data of the motor drive current;

[0093] Step S14: Based on the characteristic data of the motor drive current, perform triaxial monitoring of the vibration acceleration of the cutting head to obtain the dynamic response data of the cutting head;

[0094] Step S15: Based on the dynamic response data of the cutting head, monitor the environmental parameters of the closed-loop CNC system in laser cutting in real time to obtain environmental impact factor data;

[0095] Step S16: Based on the environmental impact factor data, timestamp the signals of each sensor to obtain time-stamped multi-source sensor data, where each sensor signal includes the optical encoder position signal, the servo motor three-phase current signal and the cutting head vibration acceleration signal.

[0096] Step S17: Generate the original system state dataset based on time-stamped multi-source sensor data;

[0097] Step S18: Extract disturbance features based on the original system state dataset to obtain system disturbance feature data, which includes steam plume disturbance feature data and servo motor thermal feature data.

[0098] In this embodiment, sensor network configuration software is selected, paired with a sensor hardware suite that supports mainstream industrial communication protocols. This software supports fine-grained configuration of parameters for optical encoders, servo motor current sensors, and vibration accelerometers. During configuration, the software's built-in frequency setting module is used to set the sampling frequency of the optical encoder to 10kHz. This parameter can be entered as 10000Hz in the "Sampling Frequency" option of the encoder property settings window and confirmed. For the servo motor current sensor, the corresponding servo motor current sensor device identifier is located in the same software interface, and its sampling frequency is set to 5kHz, i.e., 5000Hz. Similarly, the sampling frequency of the vibration accelerometer is set to 2kHz. After completing the sensor sampling frequency settings, the baud rate of the CAN bus communication protocol is further configured. In the communication protocol configuration submenu of the software, the CAN bus option is found, and the baud rate parameter is set to 1Mbps, i.e., the 1000000bps option is selected in the baud rate drop-down menu. Through the above operations, the sensor network configuration parameters are determined. For step S12, a Renishaw optical encoder is used as an example, paired with its dedicated signal acquisition and processing software. In the software, the configured optical encoder is connected and identified via the device manager. When the data acquisition function is started, the software acquires the orthogonal encoded pulse signal output by the encoder in real time at the 10kHz sampling frequency set in step S11. To improve the resolution to 0.1μm, the software's built-in 4x frequency multiplication technology module is enabled; this can be achieved simply by selecting the "4x frequency multiplication" option in the signal processing options. To address the 50Hz power frequency interference issue, the software is equipped with a hardware filtering function. By selecting the "50Hz notch filter" in the filtering settings module, power frequency interference components can be effectively filtered out, ultimately obtaining a pure original position feedback signal. In step S13, a three-phase current monitor is used, along with its dedicated current monitoring and analysis software. The monitor is connected to the three-phase power supply line of the servo motor, ensuring accurate coupling with the motor drive current signal. After the software starts, it monitors the U, V, and W three-phase currents in real time according to the configuration parameters in step S11. The three-phase current waveforms can be viewed intuitively on the software's data display interface, and the three-phase currents are converted into two-phase current components in the α-β coordinate system using the Clarke transform function module. During the transformation process, users simply need to click the "Clarke Transformation" button on the software control panel to automatically complete the coordinate transformation and obtain the motor drive current characteristic data. A MEMS triaxial accelerometer is used, along with its accompanying vibration monitoring software. The sensor is securely mounted on the X, Y, and Z axes of the cutting head to ensure the measurement range covers a dynamic response range of ±16g. In the software, the sensor's range is calibrated, selecting the ±16g range option to match the maximum acceleration changes the cutting head will encounter.The built-in signal conditioning module in the software amplifies the sensor output signal by a factor of 10 and enables low-pass filtering. The filter cutoff frequency can be set empirically to 1kHz to effectively remove high-frequency noise interference. Through these operations, the vibration acceleration response data of the cutting head during dynamic cutting can be accurately acquired. Temperature and humidity sensors and barometric pressure sensors are selected and used in conjunction with general data acquisition and environmental monitoring software. These sensors are installed inside the processing chamber of the laser cutting machine to ensure real-time sensing of environmental parameter changes. In the data acquisition software, the measurement range of the temperature and humidity sensors is set to -10°C to 60°C and 10% to 90%RH, with a sampling interval of 1 second; similarly, the sampling interval of the barometric pressure sensor is set to 1 second. After the software starts, it acquires ambient temperature, humidity, and barometric pressure data in real time, which will be used as environmental impact factor data. An industrial-grade GPS clock synchronization module, along with its clock synchronization software, is used. The GPS clock synchronization module is connected to the sensor network system. After receiving satellite signals, the module provides a microsecond-level time reference for the entire sensor network. In the software, the signals acquired by each sensor (including the optical encoder position signal, the servo motor three-phase current signal, and the cutting head vibration acceleration signal) are timestamped. Specifically, in the software's synchronization settings interface, the "GPS clock synchronization" option is selected, and the timestamp accuracy is set to 1μs. After clicking the "Start Synchronization" button, the software automatically adds a 64-bit timestamp identifier to each sensor data packet, ensuring accurate time alignment of sensor data with different sampling frequencies. The National Instruments (NI) data acquisition and processing platform can be used in conjunction with LabVIEW software. At the hardware level, NI's data acquisition card is used to aggregate the timestamped multi-source sensor data. In the LabVIEW software, the data processing flow is constructed as follows: the timestamp-marked data from each sensor is obtained from the data acquisition card through the data reading module; the data is encapsulated using the data frame structure definition module according to the format of a 2-byte header, a 1-byte data type, an 8-byte timestamp, a 2-byte data length, an N-byte payload, and a 2-byte checksum. The encapsulated standardized sensor data packets are stored in a data buffer. Then, a data fusion preprocessing module uses weighted least squares to fuse multi-sensor data from the same time point, with the weighting coefficients dynamically adjusted based on the signal-to-noise ratio of each sensor. Finally, an optimal estimation is performed using a Kalman filter to generate the original system state dataset. For a detailed implementation of step S18, please refer to the sub-steps of step S18.

[0099] Preferably, step S18 includes the following steps:

[0100] Step S181: Perform time series interpolation on the original system state dataset to obtain a system state dataset with a uniform sampling rate;

[0101] Step S182: Perform time window segmentation based on the unified sampling rate system state dataset to obtain synchronized system state data;

[0102] Step S183: Initialize the fiber optic sensor with wavelength demodulation based on the synchronization system status data to obtain the reference wavelength of the fiber optic sensor; perform spectral feature detection on the metal vapor plume based on the reference wavelength of the fiber optic sensor to obtain the plume spectral absorption data;

[0103] Step S184: Perform vapor density inversion based on plume spectral absorption data to obtain vapor plume interference characteristic data;

[0104] Step S185: Reconstruct the temperature field of the thermocouple array based on the synchronous system state data to obtain the servo motor temperature field distribution data;

[0105] Step S186: Calculate the temperature gradient of the motor windings based on the servo motor temperature field distribution data to obtain the servo motor thermal characteristic data.

[0106] In this embodiment, the `interp1` function in MATLAB is used for time series interpolation. It is assumed that the original system state dataset contains sensor data with different sampling frequencies, such as an optical encoder position signal of 10kHz, a servo motor current signal of 5kHz, and a vibration acceleration signal of 2kHz. To unify the sampling rate, 10kHz is chosen as the target sampling rate. In MATLAB, the time vector and sampled values ​​of the original data are defined; for example, for the servo motor current signal, the original time vector is `time_motor` and the sampled values ​​are `value_motor`. A unified time vector `time_unified` is created with 0.1ms intervals (the time interval corresponding to 10kHz). Cubic spline interpolation is performed using `interp1(time_motor, value_motor, time_unified, 'spline')` to obtain the interpolated current signal `value_motor_interp`. After performing the same operation on other sensor data, all interpolated data are combined into a unified sampling rate system state dataset. MATLAB is then used for time window segmentation. The time window length is set to 100ms, and the window overlap rate is 50%. The number of sampling points in each window is calculated. For a sampling rate of 10kHz, 100ms corresponds to 1000 sampling points. Data is extracted using the buffer function or by manually sliding the window. For example, for a unified sampling rate system state dataset `data`, with a window step size of 500 sampling points (corresponding to a 50% overlap rate), the segmented data matrix can be obtained using `windowedData=buffer(data,1000,500)`. To reduce spectral leakage, the Hanning window function is applied. The window function is generated using `window=hanning(1000)`, and windowing is applied to each window of data: `windowedData=windowedData×repmat(window,1,size(windowedData,2))`, finally obtaining the synchronization system state data. A tunable laser and a spectrometer are used for wavelength demodulation initialization of the fiber optic sensor. Taking the tunable laser as an example, its wavelength range covers 1530nm to 1570nm. Using the accompanying control software, the laser is set to perform wavelength scanning in steps of 0.01nm. In the software, the starting wavelength was set to 1530 nm, the ending wavelength to 1570 nm, and the scan step size to 0.01 nm. Simultaneously, a spectrometer (such as Anritsu's MS9710C) was used to monitor the spectral characteristics of the fiber optic sensor in real time. During initialization, the optical power values ​​at different wavelengths were recorded, and a wavelength-power correspondence lookup table was established. Specifically, the measurement range and resolution were set in the spectrometer software, and then the scan was started. The software automatically recorded the data and generated a lookup table to obtain the reference wavelength data for the fiber optic sensor. The Beer-Lambert law was used to invert the vapor density.The plume spectral absorption data, including light intensity attenuation information at specific wavelengths, has been obtained through step S183. Data processing is performed using OriginLab software. The light intensity attenuation data (e.g., incident light intensity I0 and transmitted light intensity I) are input, and then a nonlinear fitting is performed using the formula I=I0×exp(-α×c×l), where α is the absorption coefficient, c is the vapor concentration (density), and l is the optical path length. In OriginLab's nonlinear fitting tool, the model formula is set, and an initial guess value (e.g., α=0.1cm) is input. -1 (l=10cm), run the fitting process. The fitting result will give the value of vapor concentration c, that is, vapor density distribution data. By analyzing the spectral absorption data at different locations, the vapor plume interference characteristic data of the entire cutting area can be obtained. ANSYS Workbench is used to reconstruct the servo motor temperature field, and the temperature data measured by the thermocouple array is imported into the mechanical module of ANSYS Workbench. The geometric model and material properties (such as thermal conductivity and specific heat capacity) of the motor shell are defined using APDL (ANSYS Parametric Design Language) script. Thermocouple-measured temperature boundary conditions are applied to the model surface. Using ANSYS's finite element interpolation function, the interpolation algorithm is set to radial basis function (RBF) interpolation to ensure a smooth transition of the temperature field. After running the analysis, the software generates a three-dimensional temperature field distribution cloud map of the servo motor, providing detailed temperature distribution data, including spatial temperature gradient and time variation characteristics. The post-processing function of ANSYS Workbench is used to calculate the motor winding temperature gradient of the servo motor temperature field distribution data. In the post-processing module of ANSYS Workbench, the "Temperature Gradient" calculation tool is selected, and the calculation area is set to the motor winding part. The software automatically calculates the spatial derivative of the temperature field to obtain the temperature gradient distribution. It can extract temperature gradient values ​​at key locations (such as the contact point between the winding and the core, the winding surface, etc.) and generate a temperature gradient vector map. Simultaneously, it calculates the rate of temperature rise and, by comparing temperature field data at different times, uses ANSYS's transient analysis function to obtain the time derivative, thereby acquiring the motor's thermal characteristic data, including the rate of temperature rise and the location of the highest temperature point.

[0107] Preferably, step S2, which involves extracting interference features from the steam plume interference data, includes:

[0108] Frequency domain transformation preprocessing is performed based on steam plume interference characteristic data to obtain plume frequency domain transformation preprocessing data. The frequency domain transformation preprocessing includes windowing the steam plume interference characteristic data using a Hamming window function, with the window length set to 1024 sampling points and an overlap rate of 75%. The data length is extended to 2048 points using zero-filling technology.

[0109] Fast Fourier Transform is performed on the preprocessed data of the plume frequency domain transform to obtain the frequency domain data of the position feedback signal;

[0110] Amplitude and phase spectra are separated based on the frequency domain data of the position feedback signal to obtain separated frequency domain feature data;

[0111] Power spectral density is estimated based on the separated frequency domain characteristic data to obtain the power spectral data of the position feedback signal;

[0112] Based on the power spectrum data of the location feedback signal, spectral peak detection and location are performed to obtain frequency domain characteristic data of plume interference.

[0113] In this embodiment, MATLAB software is used for frequency domain transformation preprocessing. The steam plume interference characteristic data is imported into the MATLAB workspace. Assuming the data length is N sampling points, the Hamming window function is selected to window the data. In MATLAB, the Hamming function is used to generate a Hamming window with a length of 1024 sampling points, for example, window=hamming(1024). The data is divided into segments of 1024 sampling points each, with 75% overlap between each segment. To improve frequency resolution, zero-padding is used to extend each segment to 2048 sampling points. Specifically, 1024 zeros are added to the end of each segment, and the FFT function is used for preprocessing before the Fast Fourier Transform, ultimately obtaining the plume frequency domain transformation preprocessed data. MATLAB is then used to perform the Fast Fourier Transform (FFT). The windowed and zero-padding data obtained in step S21 is used as input, and the FFT function in MATLAB is called to perform the Fast Fourier Transform. For example, for the preprocessed data segment `data_windowed`, execute `fft_data=fft(data_windowed,2048)`, where 2048 represents the number of transformed points. The resulting `fft_data` is a complex array containing the frequency domain information of the position feedback signal. By calculating the modulus and phase angle of the complex numbers, the amplitude spectrum and phase spectrum can be obtained, specifically by `amplitude_spectrum=abs(fft_data)` and `phase_spectrum=angle(fft_data)`. MATLAB is used to separate the amplitude and phase spectra of the frequency domain data obtained in step S22. The amplitude and phase spectra are extracted from `fft_data`. The amplitude spectrum can be obtained by calculating the modulus of the complex numbers, i.e., `amplitude_spectrum=abs(fft_data)`; the phase spectrum can be obtained by calculating the phase angle of the complex numbers, i.e., `phase_spectrum=angle(fft_data)`. The amplitude and phase spectra are stored as separate arrays. For example, the amplitude spectrum and phase spectrum are stored as arrays `amplitude_spectrum` and `phase_spectrum`, respectively, using `amplitude_spectrum=abs(fft_data)` and `phase_spectrum=angle(fft_data)`. Power spectral density (PSD) estimation is then performed using MATLAB. Based on the amplitude spectrum data obtained in step S23, the Welch method is used for PSD estimation. In MATLAB, this process can be implemented using the `pwelch` function. Specifically, the window length is set to 1024 sampling points, the overlap rate is 75%, and the Hanning window function is used.For example, `window=hann(1024)`, `noverlap=768` (corresponding to a 75% overlap rate), and `nfft=2048` (FFT points). Call the `pwelch` function: `[pxx,f]=pwelch(data,window,noverlap,nfft,fs)`, where `fs` is the sampling frequency (e.g., 10kHz). The resulting `pxx` is the power spectral density estimation result, and `f` is the corresponding frequency vector. Through the above operations, the power spectral data of the position feedback signal can be obtained. Use MATLAB for spectral peak detection and location. Based on the power spectral density data `pxx` and the frequency vector `f`, set the threshold for spectral peak detection. Typically, three times the average power spectral density can be chosen as the threshold, for example, `threshold=3×mean(pxx)`. Use MATLAB's `findpeaks` function to detect spectral peaks exceeding this threshold. The specific operation is `[peaks,locs]=findpeaks(pxx,'MinPeakHeight',threshold)`, where `peaks` is the amplitude of the detected spectral peak, and `locs` is the index of the corresponding frequency point. The frequency values ​​corresponding to the spectral peaks can be obtained using f(locs). Characteristic parameters such as frequency, amplitude, and bandwidth of these spectral peaks are extracted to form frequency domain characteristic data of plume interference.

[0114] Preferably, step S2, which involves analyzing the motor parameter drift of the servo motor thermal characteristic data, includes:

[0115] Wavelet basis function selection is performed based on servo motor thermal characteristic data to obtain wavelet transform parameters. The wavelet basis function selection includes selecting Daubechies4 wavelet as mother wavelet function, setting the number of decomposition layers to 6, the scale parameter a varying from 1 to 64, and the step size of the translation parameter b set to 1 / 10 of the sampling interval.

[0116] Continuous wavelet transform is performed on the thermal characteristic data of the servo motor according to the wavelet transform parameters to obtain the time-frequency domain data of the motor temperature.

[0117] Electromagnetic constant attenuation is identified based on motor temperature time-frequency domain data to obtain raw data of electromagnetic constant attenuation.

[0118] Motor thermal characteristics are identified based on the original data of electromagnetic constant decay to obtain instantaneous parameters of motor thermal characteristics; Huang transform decomposition is performed on the instantaneous parameters of motor thermal characteristics to obtain the intrinsic mode function data of thermal characteristics.

[0119] Instantaneous frequency calculations are performed based on thermal characteristic intrinsic mode function data to obtain motor parameter drift characteristic data.

[0120] In this embodiment, MATLAB's Wavelet Toolbox is used for wavelet basis function selection and continuous wavelet transform. The servo motor thermal feature data is imported into the MATLAB workspace. The Daubechies4 wavelet is selected as the mother wavelet function, and the decomposition level is set to 6 levels. Specifically, in MATLAB, wname='db4' is used to specify the Daubechies4 wavelet, and nlevels=6 is used to set the decomposition level. The scale parameter a is set to vary from 1 to 64, and the translation parameter b has a step size of 1 / 10 of the sampling interval. For example, assuming a sampling interval of 0.1ms, the translation parameter step size is 0.01ms. The cwt function is used for continuous wavelet transform: cwt_data=cwt(thermal_data,1:64,'db4','SamplingPeriod',0.0001), where thermal_data is the thermal feature data, and 0.0001 is the sampling period (in seconds). The resulting cwt_data is the time-frequency domain data of the motor temperature. MATLAB is used to identify the electromagnetic constant decay. Based on the time-frequency domain data of motor temperature, `cwt_data`, the electromagnetic constant decay characteristics are identified by analyzing the wavelet coefficient variation trends at different scales. The specific operations are as follows: The mean and variance of the wavelet coefficients at each scale are calculated to determine the scale ranges of significant changes. For example, `mean_coeff=mean(abs(cwt_data),2)` is used to calculate the mean at each scale, and `var_coeff=var(abs(cwt_data),0,2)` is used to calculate the variance. By setting a threshold (e.g., twice the standard deviation of the mean), the scales of significant changes are identified. The frequency ranges corresponding to these scales are related to the electromagnetic constant decay. These scales and their corresponding frequency ranges are recorded to form the raw data of electromagnetic constant decay. Hilbert transform and Huang transform decomposition are performed using MATLAB. Based on the raw data of electromagnetic constant decay, a Hilbert transform is performed on the data to extract instantaneous parameters. In MATLAB, the `hilbert` function is used to transform the signal, for example, `analytic_signal=hilbert(thermal_data)`, to obtain the analytic signal. The real part of the analytic signal is the original signal, and the imaginary part is the result of the Hilbert transform. By calculating the magnitude and amplitude of the analytic signal, the instantaneous amplitude and instantaneous phase can be obtained, for example, instantaneous_amplitude=abs(analytic_signal) and instantaneous_phase=angle(analytic_signal).The instantaneous frequency is obtained by differentiating the instantaneous phase, as shown in the formula: instantaneous_frequency = diff(instantaneous_phase) / (2 × pi) × fs, where fs is the sampling frequency. Empirical Mode Decomposition (EMD) is performed, using the emd function to decompose the signal into multiple Instantaneous Mode Functions (IMFs). For example, imf = emd(instantaneous_frequency), resulting in multiple IMF components. Each IMF component represents the characteristics of a different frequency band. MATLAB is used to calculate the instantaneous frequency of the thermal characteristic IMF data. Based on the IMF component imf obtained in step S28, a Hilbert transform is performed on each IMF component to obtain the instantaneous frequency of each IMF. Specifically, for each IMF component imf_i, the hilbert function is used to obtain the analytic signal analytic_signal_i = hilbert(imf_i). The instantaneous frequency is then calculated.

[0121] `instantaneous_frequency_i = diff(angle(analytic_signal_i)) / (2 × pi) × fs`, where `fs` is the sampling frequency. Extract the instantaneous frequency of each IMF component and calculate its characteristic parameters such as average frequency, frequency variance, and frequency drift rate. For example:

[0122] mean_freq_i = mean(instantaneous_frequency_i), var_freq_i = var(instantaneous_frequency_i), drift_rate_i = diff(instantaneous_frequency_i) / dt, where dt is the time interval. These characteristic parameters will serve as motor parameter drift characteristic data.

[0123] Of particular importance, the identification of electromagnetic constant decay based on motor temperature time-frequency domain data also includes the following steps:

[0124] Wavelet energy spectral density was calculated from the time-frequency domain data of motor temperature to obtain wavelet energy spectral distribution data;

[0125] The energy of the preset sensitive frequency band is integrated based on the wavelet energy spectrum distribution data to obtain the frequency band energy integral data.

[0126] Time series trend identification is performed based on frequency band energy integral data to obtain temperature change trend characteristic data;

[0127] Based on the temperature change trend characteristic data, the attenuation coefficient is fitted with the preset temperature-electromagnetic constant attenuation mapping model to obtain preliminary attenuation coefficient data;

[0128] The initial attenuation coefficient data is subjected to sliding window mean filtering to obtain smooth attenuation coefficient data;

[0129] The cumulative amount is calculated based on the smooth attenuation coefficient data to obtain the cumulative electromagnetic constant attenuation data; the data format is standardized based on the cumulative electromagnetic constant attenuation data to obtain the original electromagnetic constant attenuation data.

[0130] In this embodiment, the Wavelet Toolbox in MATLAB is used to calculate the wavelet energy spectral density. Based on the time-frequency domain data of the motor temperature, cwt_data, the absolute values ​​of the wavelet coefficients are calculated using the abs function, and then squared to obtain the wavelet energy spectral density. Specifically, wavelet_energy = abs(cwt_data).^2. After generating the wavelet energy spectral distribution data, the energy distribution at different frequencies and time points can be visualized by plotting two-dimensional or three-dimensional heatmaps. This visualization effect is achieved using the imagesc or surf functions in MATLAB. Bandwidth energy integration is performed using MATLAB. Based on the wavelet energy spectral distribution data wavelet_energy, a sensitive frequency band range is preset, for example, selecting a band related to the temperature change of the motor windings as 1-50Hz. The integration range is defined using freq_vector (frequency vector) and time_vector (time vector), and the wavelet energy spectrum is integrated within this frequency band. In MATLAB, the cumtrapz function can be used for numerical integration: energy_integral = cumtrapz(freq_vector, wavelet_energy, 1), where dimension 1 corresponds to the frequency axis. The obtained `energy_integral` is the bandgap energy integral data. MATLAB's Signal Processing Toolbox is used for time series trend identification. Based on the obtained bandgap energy integral data `energy_integral`, the data is detrended to eliminate linear or nonlinear trends. The `detrend` function is used to remove linear trends, for example:

[0131] detrended_data=detrend(energy_integral). Apply a moving average filter to smooth the data using the movmean function, for example, smoothed_data=movmean(detrended_data,[5 5]), the window size can be adjusted according to the data characteristics. Perform time series decomposition to identify trend, seasonality and residual components. Use the stl function (seasonal decomposition) to decompose the data, decomposed_data=stl(smoothed_data,'Seasonality',10), where the seasonal period is set to 10. By analyzing the trend components after decomposition, temperature change trend feature data is obtained. Use MATLAB's Curve Fitting Toolbox to fit the attenuation coefficient. Based on the temperature change trend feature data, combine the preset temperature-electromagnetic constant attenuation mapping model. Assuming the mapping model is a linear relationship, the fit function can be used for linear fitting. For example, fitmodel=fit(time_vector',trend_data','poly1'), where time_vector is the time vector, trend_data is the temperature change trend data, and poly1 represents a first-order polynomial fitting. The obtained fitting model fitmodel contains slope and intercept parameters, and the slope is the initial decay coefficient data. If the model is nonlinear, the appropriate fitting type can be selected according to the actual relationship, such as the exponential decay model fitmodel=fit(time_vector',trend_data','exp1'). MATLAB is used for sliding window mean filtering. Based on the initial decay coefficient data initial_decay_coeff, a suitable sliding window size is selected, for example, the window size is set to 10 data points for smoothing. In MATLAB, the movmean function is used to perform sliding window mean filtering on the data, smoothed_decay_coeff=movmean(initial_decay_coeff,

[55] ), and the window parameter

[55] indicates that 5 points are taken before and after the current point as the center for averaging. The obtained smoothed_decay_coeff is the smoothed decay coefficient data. MATLAB is used for cumulative calculation and data format standardization. Based on the smoothed decay coefficient data smoothed_decay_coeff, the cumulative decay is calculated. The cumsum function can be used to calculate the cumulative sum: cumulative_decay = cumsum(smoothed_decay_coeff), which yields the cumulative data of electromagnetic constant decay.To standardize the data format, the cumulative decay data is normalized to a preset range, such as 0-1. The `rescale` function, `standardized_decay_data=rescale(cumulative_decay,0,1)`, is used to obtain the original data on electromagnetic constant decay.

[0132] Preferably, the construction of the system parameter drift indicator library based on the plume interference frequency domain indicator data and the motor parameter drift indicator data in step S2 includes:

[0133] Construct a plume interference feature vector based on plume interference frequency domain characteristic data;

[0134] Construct a motor drift feature vector based on motor parameter drift symptom data;

[0135] Multidimensional feature space mapping is performed based on the plume interference feature vector and the motor drift feature vector to obtain the fused system state feature space data;

[0136] Clustering and pattern recognition are performed on the fusion system state feature space data to obtain a system parameter drift sign library.

[0137] In this embodiment, MATLAB is used to construct the plume interference feature vector. The plume interference frequency domain characteristic data obtained from step S25 includes information such as spectral peak frequency, amplitude, and bandwidth. For example, suppose four significant spectral peaks are detected, each with three feature parameters: frequency, amplitude, and bandwidth. In MATLAB, these feature parameters can be combined into a 12-dimensional feature vector. The specific operation is as follows: create a 1×12 array, where the first four elements are the frequencies of each spectral peak, the next four elements are the corresponding amplitudes, and the last four elements are the bandwidths. For example:

[0138] `vapor_vector=[f1,f2,f3,f4,a1,a2,a3,a4,b1,b2,b3,b4]`, where f, a, and b represent frequency, amplitude, and bandwidth, respectively. A motor drift feature vector is constructed using MATLAB. The motor parameter drift characteristic data obtained in step S29 includes the mean, variance, and drift rate of the instantaneous frequency; the mean, variance, and drift rate of the amplitude; and the mean and variance of the phase, totaling eight feature parameters. In MATLAB, these parameters are combined into an 8-dimensional feature vector. The specific operation is as follows: Create a 1×8 array to store the mean frequency, frequency variance, frequency drift rate, mean amplitude, amplitude variance, amplitude drift rate, mean phase, and phase variance in sequence. For example;

[0139] `motor_vector = [f_mean, f_var, f_drift, a_mean, a_var, a_drift, phase_mean, phase_var]`. Multidimensional feature space mapping is performed using Python and its NumPy and scikit-learn libraries. The plume disturbance feature vector `vapor_vector` and the motor drift feature vector `motor_vector` obtained in step S211 are fused. The two vectors are concatenated into a 20-dimensional original feature vector. For example, `combined_vector = np.concatenate((vapor_vector, motor_vector))`. Principal component analysis (PCA) is used for dimensionality reduction, retaining 95% of the variance contribution rate. In Python, the PCA class is imported and initialized: `pca = PCA(n_components = 0.95)`. The concatenated feature vector is fitted and transformed: `reduced_vector = pca.fit_transform(combined_vector.reshape(1, -1))`. The resulting `reduced_vector` is the dimensionality-reduced fused system state feature space data. Clustering and pattern recognition are performed using Python and its scikit-learn library. Based on fused system state feature space data, K-means clustering analysis is used. The number of cluster centers is determined, for example, 8. In Python, the KMeans class is imported and initialized: `kmeans = KMeans(n_clusters = 8)`. Clustering training is performed on the fused feature data: `kmeans.fit(fused_features)`, where `fused_features` is the fused feature dataset. After training, each data point is assigned to a cluster center. By analyzing the clustering results, a system parameter drift signature database is established. Each cluster center represents a typical parameter drift pattern, and the signature database will be used for subsequent real-time monitoring and pattern matching.

[0140] Preferably, step S3 includes the following steps:

[0141] Step S31: Initialize the system parameters using the recursive least squares method based on the system parameter drift symptom library to obtain the initial parameters of the recursive algorithm;

[0142] Step S32: Real-time data acquisition is performed on the closed-loop CNC system in laser cutting to obtain system input and output data; regression matrix is ​​constructed on the system input and output data based on the initial parameters of the recursive algorithm to obtain the regression matrix of the closed-loop CNC system;

[0143] Step S33: Perform recursive least squares parameter updates based on the regression matrix of the closed-loop CNC system to obtain real-time CNC system parameter estimation data;

[0144] Step S34: Extract the feedforward compensation coefficients based on the parameter estimation data of the real-time CNC system to obtain the original data of the feedforward compensation coefficients;

[0145] Step S35: Calculate the forgetting factor based on the original data of the feedforward compensation coefficient to obtain the forgetting factor data of the CNC system;

[0146] Step S36: Adjust the weights of the original data of the feedforward compensation coefficients based on the forgetting factor data of the CNC system to obtain the weighted feedforward compensation data of the CNC system;

[0147] Step S37: Perform parameter oscillation suppression on the feedforward compensation data of the weighted CNC system to obtain stable feedforward compensation parameters.

[0148] In this embodiment, the recursive least squares method initialization is performed using MATLAB's System Identification Toolbox. The dimension of the parameter vector is set to 6 via the graphical interface or command-line window, and the parameter estimate θ(0) is initialized as a 6-dimensional column vector [1;0;0;1;0;0]. Simultaneously, the covariance matrix P(0) is initialized to a 6×6 identity matrix multiplied by 100. The forgetting factor λ is set to 0.95. These initial parameter settings form the basis for subsequent recursive least squares calculations. Real-time data from the closed-loop CNC system in laser cutting is collected from the sensor network using MATLAB's data import function. A regression matrix is ​​constructed by reading the system input / output data stored in the workspace. In MATLAB, combining the initial parameters of the recursive algorithm obtained in step S31, the output position signal and input control signal are sequentially arranged and combined into a regression matrix φ(k), which contains the values ​​of y(k-1), y(k-2), u(k-1), u(k-2), u(k-3), and u(k-4). This matrix will be used for subsequent recursive least squares parameter updates. Parameter updates are performed using MATLAB's recursive least squares method. The regression matrix and real-time system input / output data are input for recursive calculation. MATLAB automatically iterates and updates the gain matrix, parameter estimates, and covariance matrix according to the recursive least squares formula. This yields real-time CNC system parameter estimates. Feedforward compensation coefficients are extracted from the real-time CNC system parameter estimates using MATLAB's vector manipulation functions. Based on the previously obtained parameter estimate vector θ, the 3rd to 6th elements related to the input signal are extracted as feedforward compensation coefficients. These coefficients are saved as the original feedforward compensation coefficient data. The forgetting factor is calculated using MATLAB. Based on the original feedforward compensation coefficient data K_ff, the dynamic forgetting factor is calculated using MATLAB's mathematical functions. A variable forgetting factor formula is used, combined with the set minimum and maximum forgetting factors and a time constant, to calculate the forgetting factor values ​​for different iteration numbers. MATLAB is used to adjust the weights of the original feedforward compensation coefficient data. The original feedforward compensation coefficient data K_ff is combined with the dynamic forgetting factor and processed using MATLAB's weighted average function. Using the dynamic forgetting factor obtained in the previous steps and the weighted feedforward compensation data from the previous time step, the weighted feedforward compensation data for the current time step of the CNC system is calculated. This step is implemented using MATLAB's vector operation function to ensure that the weighted feedforward compensation data can be effectively used for subsequent parameter oscillation suppression and system stability improvement. For detailed implementation of step S37, please refer to the sub-steps of step S37.

[0149] Preferably, step S37 includes the following steps:

[0150] Step S371: Extract historical data from the sliding window based on the weighted CNC system feedforward compensation data to obtain a set of historical CNC system data windows;

[0151] Step S372: Perform weighted least squares fitting based on the historical CNC system data window set to obtain the feedforward compensation coefficients of the CNC system;

[0152] Step S373: Calculate the parameter change rate based on the feedforward compensation coefficient of the CNC system to obtain the compensation coefficient change rate data;

[0153] Step S374: Perform threshold detection and amplitude limiting based on the rate of change of the compensation coefficient to obtain the corrected feedforward compensation coefficient;

[0154] Step S375: Perform system stability assessment based on the modified feedforward compensation coefficient to obtain system stability assessment data; select an adjustment strategy based on the system stability assessment data and the preset adjustment strategy set to obtain the effective strategy adjustment parameters of the system;

[0155] Step S376: Adjust the parameters based on the effective strategy of the system to perform parameter convergence verification and obtain the convergence verification data of the closed-loop CNC system;

[0156] Step S377: Confirm the final parameters based on the convergence verification data of the closed-loop CNC system to obtain stable feedforward compensation parameters.

[0157] In this embodiment, MATLAB's sliding window function is used for historical data extraction. Based on the weighted CNC system feedforward compensation data, the sliding window length is set to 50 sampling points to capture recent system behavior. In MATLAB, the feedforward compensation coefficients for the past 50 moments are extracted using the sliding window operation, forming a 50×4 historical data matrix. For example, the data can be extracted using the `buffer` function or by manually sliding the window to obtain a set of historical CNC system data windows. MATLAB is used for weighted least squares fitting. Based on the historical data window set, linear fitting is performed on the data within each window. In MATLAB, the weight coefficients are set to an exponentially decreasing form, for example, weight coefficient w(i) = exp(-(50-i) / 10), where i is the position of the data point in the window. The trend of the feedforward compensation coefficients is fitted using the weighted least squares method to obtain the CNC system feedforward compensation coefficients. MATLAB is used to calculate the rate of change of parameters. Based on the CNC system feedforward compensation coefficients, the rate of change of the compensation coefficients at adjacent moments is calculated. In MATLAB, the rate of change is calculated using difference calculations, for example, dK_ff = diff(K_ff) calculates the first-order difference to obtain the rate of change data of the compensation coefficient. Simultaneously, the second derivative is calculated to assess the acceleration of the change, for example, d2K_ff = diff(K_ff, 2). MATLAB is then used for threshold detection and limiting. Based on the rate of change data of the compensation coefficient, a threshold of ±0.1 / s is set. In MATLAB, by comparing the rate of change with the threshold, data points exceeding the threshold are identified. For cases exceeding the threshold, limiting processing is performed, for example, restricting the portion exceeding the threshold to within the threshold range. Specifically, if the rate of change exceeds 0.1 / s, the corrected coefficient is the coefficient of the previous time step plus 0.1 × Ts (Ts is the sampling period); if it is below -0.1 / s, 0.1 × Ts is subtracted. The corrected feedforward compensation coefficient is obtained. MATLAB is then used for system stability evaluation. Based on the corrected feedforward compensation coefficient, the statistical characteristics of the compensation coefficient are calculated, including mean, variance, skewness, and kurtosis. In MATLAB, statistical analysis tools are used to calculate these indicators; for example, `mean(K_ff_corrected)` calculates the mean, and `var(K_ff_corrected)` calculates the variance. The statistical stability of the coefficient sequence is assessed through a normality test. Based on the system stability assessment data and a pre-defined set of adjustment strategies (such as conservative, standard, and aggressive strategies), an appropriate adjustment strategy is selected. For example, a conservative strategy is chosen when the variance is less than 0.01; a standard strategy is chosen when the variance is between 0.01 and 0.1; and an aggressive strategy is chosen when the variance is greater than 0.1. The system's effective strategy adjustment parameters are obtained. MATLAB is used to test the parameter convergence. Based on the system's effective strategy adjustment parameters, the Lyapunov exponent of the parameter sequence is calculated to assess convergence.In MATLAB, the Lyapunov exponent calculation function, such as `lyapunov_exponent`, is used in the Nonlinear Time Series Analysis Toolbox. Test parameters, such as embedding dimension and delay time, are set to obtain the Lyapunov exponent. A negative exponent indicates convergence of the parameter sequence. Simultaneously, the 95% confidence interval of the parameter estimate is calculated to verify the estimation accuracy. Through these operations, convergence verification data for the closed-loop CNC system is obtained. Final parameter confirmation is then performed using MATLAB. Based on the convergence verification data, it is checked whether the parameter change rate is less than 0.01 and the Lyapunov exponent is less than -0.1 over 10 consecutive sampling periods. In MATLAB, this verification process is implemented through logical judgment and conditional checks. If the conditions are met, parameter convergence is confirmed, and the final feedforward compensation parameter vector is output. This step is implemented using MATLAB's conditional statements and logical operations to ensure that final confirmation is only performed when the parameters are stable and convergent, thus obtaining stable and reliable feedforward compensation parameters.

[0158] Preferably, step S4 includes the following steps:

[0159] Step S41: Construct a mathematical model of the closed-loop system based on the stable feedforward compensation parameters to obtain the system state-space model;

[0160] Step S42: Extract the characteristic polynomial coefficients based on the system state-space model to obtain the characteristic polynomial parameters of the closed-loop system;

[0161] Step S43: Identify the system pole locations based on the characteristic polynomial parameters of the closed-loop system to obtain the closed-loop system pole data;

[0162] Step S44: Identify stability criteria based on the closed-loop system pole data to obtain the stability coefficients of the closed-loop system;

[0163] Step S45: Evaluate the stability margin of the closed-loop CNC system in laser cutting based on the stability coefficient of the closed-loop system to obtain the stability margin of the closed-loop CNC system.

[0164] Step S46: Iteratively optimize the speed loop gain parameters based on the stability margin of the closed-loop CNC system to obtain the optimized speed loop control parameters.

[0165] In this embodiment, the Control System Toolbox of MATLAB is used to construct the mathematical model of the closed-loop system. Based on the stable feedforward compensation parameters, the system's state matrix, input matrix, output matrix, and direct transfer matrix are defined. For example, the state matrix A is a 4×4 matrix, the input matrix B is a 4×1 matrix, the output matrix C is a 1×4 matrix, and the direct transfer matrix D is 0. In MATLAB, the ss function is used to create the state-space model: sys=ss(A,B,C,D). The characteristic polynomial coefficients are extracted using the Control System Toolbox of MATLAB. Based on the system state-space model sys, the tf function is used to convert it into transfer function form: [num,den]=tf(sys). The resulting denominator polynomial den is the characteristic polynomial parameter of the closed-loop system. The pole location of the system is identified using the pole function of MATLAB. Based on the characteristic polynomial parameter den, the pole location of the system is directly calculated using the pole(sys) function. For example, poles=pole(sys) will return a complex vector containing the pole locations of the system. The stability criteria are identified using the tf and tf2ss functions of MATLAB. Based on the pole locations (poles), a Routh-Hurwitz table is constructed to evaluate the system's stability. The `tf` function is used to obtain the system's transfer function form, and then the `tf2ss` function is used to convert it to state-space form. The `roots` function is used to find the roots of the characteristic equation, and a Routh-Hurwitz table is constructed based on these roots. The coefficients in the Routh-Hurwitz table are the closed-loop system stability coefficients. Specifically, `roots(den)` obtains the roots of the characteristic equation. The stability of the system can be determined by analyzing the number of sign changes in the first column of the Routh table. For example, if there are no sign changes in the first column, the system is stable. The stability margin is evaluated using MATLAB's `margin` function. Based on the obtained Routh table coefficients, the system's gain margin and phase margin are calculated. In MATLAB, the `margin` function is used to plot a Bode plot and automatically calculate the gain margin and phase margin. For example, `margin(sys)` will display the Bode plot and label the specific values ​​of the gain margin and phase margin on the plot. For detailed implementation of step S46, please refer to the sub-steps of step S46.

[0166] Of particular importance, step S46 includes the following steps:

[0167] Step S461: Construct the optimization objective function based on the stability margin of the closed-loop CNC system to obtain the speed loop optimization objective function parameters; perform gradient calculation and direction determination based on the speed loop optimization objective function parameters to obtain the speed loop gain optimization gradient vector;

[0168] Step S462: Perform adaptive step size calculation based on the velocity loop gain optimization gradient vector to obtain velocity loop gain optimization step size data; update the velocity loop gain parameters based on the velocity loop gain optimization step size data to obtain candidate velocity loop gain parameters;

[0169] Step S463: Perform constraint condition verification based on the candidate parameters of the velocity loop gain to obtain velocity loop gain constraint verification data;

[0170] Step S464: Perform system response simulation based on velocity loop gain constraint verification data to obtain closed-loop system simulation response characteristic data; evaluate performance indicators based on closed-loop system simulation response characteristic data to obtain velocity loop performance evaluation data.

[0171] Step S465: Based on the velocity loop performance evaluation data, perform convergence judgment and iterative control to obtain gain optimization iterative control data;

[0172] Step S466: Confirm the optimal parameters based on the gain optimization iterative control data to obtain the optimal speed loop gain data; construct a real-time adjustment mechanism based on the optimal speed loop gain data to obtain the speed loop optimization control parameters.

[0173] In this embodiment, the Optimization Toolbox of MATLAB is used to construct the optimization objective function. Based on the stability margin of the closed-loop CNC system, the performance index function J=∫0^∞[e 2 (t)+ρ1u 2 (t)+ρ2u̇ 2The expression is: ∫[e(t)]dt, where e(t) is the position tracking error, u(t) is the control input, u̇(t) is the rate of change of the control input, and the weighting coefficients are ρ1=0.1 and ρ2=0.01. In MATLAB, the fmincon function is used to construct the optimization objective function. The objective function is defined as the weighted sum of the squares of the position tracking error, the squares of the control input, and the squares of the rate of change of the control input. By calculating the partial derivatives of the objective function with respect to the velocity loop gain parameters, the velocity loop gain optimization gradient vector ∇J=[∂J / ∂Kp,∂J / ∂Ki,∂J / ∂Kd] is obtained. The step size adaptive solution is performed using MATLAB's Optimization Toolbox. Based on the obtained gradient vector ∇J, the Armijo criterion is used to determine the optimal step size αk. In MATLAB, the fminunc function is used for unconstrained optimization, setting the initial step size to 1, and the step size satisfying the conditions is determined through backtracking search. For example, optimization options are set using `options=optimoptions('fminunc','Algorithm','quasi-newton','Display','iter')`, and optimization is performed by calling `[x,fval]=fminunc(fun,x0,options)`. The velocity loop gain parameters are updated iteratively to obtain candidate velocity loop gain parameters Kp, Ki, and Kd. MATLAB is used to verify the constraints. Based on the candidate velocity loop gain parameters, it is verified whether the parameters are within the design range. In MATLAB, the constraints are set as follows: proportional gain 0.1≤Kp≤10, integral gain 0.01≤Ki≤1, and derivative gain 0.001≤Kd≤0.1. Conditional statements (such as if statements) are used to check whether each parameter is out of range; parameters out of range are projected onto boundary values. For example, if Kp>10, Kp=10; if Ki<0.01, Ki=0.01. System response simulation is performed using MATLAB's Simulink. Based on the obtained constraint verification data, a Simulink model of the closed-loop system is constructed. In the model, a unit step signal is input, the simulation time is set to 1 second, and the sampling interval is 0.1 ms. The simulation yields response characteristic data such as the system's rise time (tr), settling time (ts), overshoot (Mp), and steady-state error (ess). In MATLAB, the simulation model is run using the `sim` function, and the simulation results are obtained using the `get_param` function. Based on these response characteristic data, the system's performance indicators are evaluated, and velocity loop performance evaluation data is obtained. MATLAB is used for convergence judgment and iterative control. Based on the performance evaluation data, the convergence condition is set as a performance score change rate of less than 0.1% over 5 consecutive iterations or an iteration count exceeding 100. Iterative control is implemented in MATLAB using loop structures (such as for loops or while loops).In each iteration, the current performance score is compared with the scores of the previous iterations, and the rate of change is calculated. If the convergence condition is met, the iteration stops; otherwise, the next iteration continues. MATLAB is used to confirm the optimal parameters and construct the real-time adjustment mechanism. Based on the obtained iterative control data, the parameter combination with the highest performance score is selected as the optimal solution from the iteration history. In MATLAB, the optimal velocity loop gain parameters are found by comparing the performance scores at different iteration numbers. The robustness of the system under this parameter combination is verified by conducting 1000 random parameter perturbation tests using the Monte Carlo method. Based on the optimal velocity loop gain data, a real-time adjustment mechanism is constructed, for example, increasing the proportional gain Kp when the system error exceeds a threshold, increasing the integral gain Ki when a steady-state error exists, and decreasing the derivative gain Kd when the system oscillates.

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

[0175] Step S51: Acquire the original trajectory data of the closed-loop CNC system in laser cutting to obtain the original cutting trajectory data; Based on the speed loop optimization control parameters, reconstruct the original cutting trajectory data to obtain the standard cutting trajectory data.

[0176] Step S52: Calculate the curvature based on the standard cutting trajectory data to obtain the cutting trajectory curvature distribution data;

[0177] Step S53: Extract geometric feature parameters based on the curvature distribution data of the cutting trajectory to obtain the geometric feature data of the cutting trajectory;

[0178] Step S54: Calculate the rate of change of tangent direction based on the geometric feature data of the cutting trajectory to obtain the feature data of the change of cutting direction;

[0179] Step S55: Perform velocity planning and prediction based on the cutting direction change feature data to obtain the cutting predicted velocity distribution data;

[0180] Step S56: Real-time trajectory data acquisition is performed on the closed-loop CNC system during laser cutting to obtain the current cutting trajectory data;

[0181] Step S57: Perform predictive compensation control on the current cutting trajectory data based on the cutting prediction speed distribution data to obtain the final cutting trajectory control data.

[0182] In this embodiment, the CNC system built into the laser cutting equipment or a separate data acquisition device is used to record the real-time position information of the cutting head during the processing with high precision. This raw trajectory data typically includes the coordinates of the cutting path and the corresponding timestamps. Based on the control parameters optimized by the speed loop, the raw trajectory data is recalculated and adjusted to ensure the cutting trajectory meets new speed and precision requirements. Curvature calculation can be performed using mathematical software such as MATLAB or the NumPy library in Python. Standard cutting trajectory data is imported; this data is typically a series of coordinate points. The curvature of the trajectory is obtained by calculating the vector change between adjacent points. Specifically, the trajectory curve is parameterized, the tangent vector at each point is calculated, and the rate of change of the tangent vector is analyzed to obtain the curvature distribution data of the trajectory. Based on the curvature distribution data of the cutting trajectory, the Shapely library in Python or the CurveFitting Toolbox in MATLAB is used to analyze the trajectory curve and extract features such as maximum curvature, average curvature, and rate of change of curvature. In step S54, when calculating the rate of change of the tangent direction, mathematical software such as MATLAB or Python can be used. Based on the curvature distribution data of the cutting trajectory, the rate of change of the tangent direction is obtained by calculating the angular change of the tangent direction between adjacent points. Specifically, at each point on the trajectory curve, the angle between the tangent direction and a reference direction (such as the X-axis) is calculated, and the rate of change of the angle between adjacent points is obtained. Based on the tangent direction change characteristic data and combined with the dynamic characteristics of the laser cutting equipment, the appropriate cutting speed for different trajectory segments is predicted. This involves a comprehensive consideration of the equipment's acceleration, deceleration, and maximum speed parameters. Through simulation and optimization rules, a speed distribution curve matching the trajectory geometry is generated, i.e., the predicted cutting speed distribution data. Real-time trajectory data acquisition can be achieved through position sensors and data acquisition cards on the laser cutting equipment. These sensors monitor the position of the cutting head in real time and transmit the data to the control system. Using data acquisition and analysis software such as MATLAB or LabVIEW, the current position information of the cutting head can be acquired and recorded in real time, forming the current cutting trajectory data. For a detailed implementation process of step S57, please refer to the sub-steps of step S57.

[0183] Of particular importance, step S57 includes the following steps:

[0184] Step S571: Calculate motion parameters based on the predicted cutting velocity distribution data to obtain system dynamics characteristic data; train a support vector regression model based on the system dynamics characteristic data to obtain a laser cutting accuracy prediction model;

[0185] Step S572: Based on the laser cutting accuracy prediction model, predict the accuracy deviation of the current cutting trajectory data to obtain the cutting accuracy deviation prediction data;

[0186] Step S573: Calculate the feedforward compensation amount based on the cutting accuracy deviation prediction data to obtain the cutting trajectory compensation vector;

[0187] Step S574: Correct the current cutting trajectory data based on the cutting trajectory compensation vector to obtain corrected trajectory coordinate data; perform a smoothness check based on the corrected trajectory coordinate data to obtain cutting trajectory smoothness verification data;

[0188] Step S575: Based on the cutting trajectory smoothness verification data, the compensation amount is limited to obtain the compensation amount constraint data;

[0189] Step S576: Evaluate the trajectory quality based on the compensation constraint data to obtain the cutting trajectory quality evaluation data; confirm the final trajectory based on the cutting trajectory quality evaluation data to obtain the final cutting trajectory control data.

[0190] In this embodiment, Python's NumPy and scikit-learn libraries are used to calculate motion parameters and train the model. Based on the predicted velocity distribution data, the system's dynamic characteristics, including acceleration and angular acceleration, are calculated. Assuming the velocity distribution data is stored as an array, acceleration is calculated using NumPy:

[0191] `acceleration = np.gradient(speed_distribution, time_interval)`, where `time_interval` is the sampling time interval. Based on the calculated dynamic feature data, the SVR (Support Vector Regression) model from scikit-learn is used for training. Specifically, the SVR class is imported, and the model is initialized as `svr_model = SVR(kernel='rbf', C=100, gamma=0.1)`. The model is trained using the dynamic feature data as input and the cutting accuracy deviation as output: `svr_model.fit(X_train, y_train)`, where `X_train` is the input feature of the training set, and `y_train` is the corresponding output label. The trained SVR model is then used to predict the accuracy deviation of the current cutting trajectory data. The current cutting trajectory data has been preprocessed and feature extracted, forming a feature matrix `X_current` in the same format as the training set. `X_current` is input into the SVR model: `predicted_deviation = svr_model.predict(X_current)`, obtaining the predicted cutting accuracy deviation data `predicted_deviation`. This predicted data contains the deviation that will occur at each trajectory point. The feedforward compensation is calculated based on the predicted deviation data. The calculation of the feedforward compensation is based on the magnitude and direction of the deviation, as well as the dynamic characteristics of the system. Assuming the compensation coefficient matrix is ​​known, matrix operations can be used to calculate the compensation. For example, in Python, matrix multiplication can be performed using NumPy: `compensation_vector = compensation_matrix @ predicted_deviation`, where `compensation_matrix` is the pre-calibrated compensation coefficient matrix, and `predicted_deviation` is the predicted deviation vector. The resulting `compensation_vector` is the cutting trajectory compensation vector, indicating the amount of compensation to be applied at each trajectory point to offset the predicted deviation. The current cutting trajectory data is then corrected based on the compensation vector. This involves superimposing the original trajectory data with the compensation vector to generate corrected trajectory coordinate data. Assuming the original trajectory coordinates are stored as a two-dimensional array `original_trajectory`, and the compensation vector is `compensation_vector`, the corrected trajectory coordinate data can be obtained through simple addition: `corrected_trajectory = original_trajectory + compensation_vector`. Smoothing functions from Python's SciPy library can be used to ensure the smoothness of the corrected trajectory.The smoothed trajectory is verified to check if it meets processing requirements, such as continuity and smoothness, resulting in cutting trajectory smoothness verification data. Compensation amount limits are then imposed on the smoothed trajectory data. The purpose of compensation amount limits is to prevent over-compensation from causing trajectory distortion. Conditional statements can be used to implement compensation amount limits. For example, in Python, upper and lower limits for compensation amount can be set, and each compensation value can be judged and adjusted: if the compensation value exceeds the upper limit, it is set to the upper limit value; if it is below the lower limit, it is set to the lower limit value. The compensation amount data processed in this way is the compensation amount constraint data. The compensation amount constraint data is used to perform quality assessment on the corrected trajectory coordinate data. Quality assessment can be based on the statistical characteristics of trajectory deviation, such as root mean square error (RMS). In Python, the RMS value can be calculated using NumPy: rms_error=np.sqrt(np.mean(compensated_trajectory×2)), where compensated_trajectory is the compensated trajectory deviation data. Based on the RMS value and other quality indicators, a comprehensive quality evaluation of the cutting trajectory is performed. If the trajectory quality meets the preset standard, such as an RMS error of less than 0.01mm, then the trajectory is confirmed as the final cutting trajectory control data.

[0192] Therefore, the embodiments should be considered as exemplary and non-limiting in all respects, and the scope of the invention is defined by the appended claims rather than the foregoing description. Thus, it is intended that all variations falling within the meaning and scope of the equivalents of the application be incorporated into the invention.

[0193] The above description is merely a specific embodiment of the present invention, enabling those skilled in the art to understand or implement the invention. Various modifications to these embodiments will be readily apparent to those skilled in the art, and the general principles defined herein may be implemented in other embodiments without departing from the spirit or scope of the invention. Therefore, the present invention is not to be limited to the embodiments shown herein, but is to be accorded the widest scope consistent with the principles and novel features of the invention herein.

Claims

1. A method for automatically detecting control model parameters of a closed-loop CNC system in laser cutting, characterized in that, Includes the following steps: Step S1: Obtain the original system state dataset; extract disturbance features based on the original system state dataset to obtain system disturbance feature data, which includes steam plume disturbance feature data and servo motor thermal feature data; Step S2: Extract interference features from the steam plume interference characteristic data to obtain plume interference frequency domain characteristic data; perform motor parameter drift analysis on the servo motor thermal characteristic data to obtain motor parameter drift characteristic data; construct a system parameter drift characteristic library based on the plume interference frequency domain characteristic data and the motor parameter drift characteristic data. Step S3: Based on the system parameter drift symptom library, perform online identification and weight adjustment of the system feedforward compensation coefficients to obtain weighted CNC system feedforward compensation data; perform parameter oscillation suppression on the weighted CNC system feedforward compensation data to obtain stable feedforward compensation parameters; Step S4: Quantify the system stability margin based on the stable feedforward compensation parameters to obtain the stability margin of the closed-loop CNC system; perform iterative optimization of the speed loop gain parameters based on the stability margin of the closed-loop CNC system to obtain the optimized speed loop control parameters. Step S5: Predict the cutting speed based on the speed loop optimization control parameters to obtain the predicted cutting speed distribution data; perform predictive compensation control on the current cutting trajectory data based on the predicted cutting speed distribution data to obtain the final cutting trajectory control data.

2. The method for automatic detection of control model parameters of a closed-loop CNC system in laser cutting according to claim 1, characterized in that, Step S1 includes the following steps: Step S11: Initialize the sensor network configuration of the closed-loop CNC system in laser cutting to obtain the sensor network configuration parameters. The sensor network initialization configuration includes setting the sampling frequency of the optical encoder to 10kHz, the sampling frequency of the servo motor current sensor to 5kHz, and the sampling frequency of the vibration accelerometer to 2kHz. At the same time, configure the CAN bus communication protocol baud rate to 1Mbps. Step S12: Acquire position signals from the optical encoder based on the sensor network configuration parameters to obtain the raw position feedback signal; Step S13: Monitor the three-phase current of the servo motor based on the original position feedback signal to obtain the characteristic data of the motor drive current; Step S14: Based on the characteristic data of the motor drive current, perform triaxial monitoring of the vibration acceleration of the cutting head to obtain the dynamic response data of the cutting head; Step S15: Based on the dynamic response data of the cutting head, monitor the environmental parameters of the closed-loop CNC system in laser cutting in real time to obtain environmental impact factor data; Step S16: Based on the environmental impact factor data, timestamp the signals of each sensor to obtain time-stamped multi-source sensor data, where each sensor signal includes the optical encoder position signal, the servo motor three-phase current signal and the cutting head vibration acceleration signal. Step S17: Generate the original system state dataset based on time-stamped multi-source sensor data; Step S18: Extract disturbance features based on the original system state dataset to obtain system disturbance feature data, which includes steam plume disturbance feature data and servo motor thermal feature data.

3. The method for automatic detection of control model parameters of a closed-loop CNC system in laser cutting according to claim 2, characterized in that, Step S18 includes the following steps: Step S181: Perform time series interpolation on the original system state dataset to obtain a system state dataset with a uniform sampling rate; Step S182: Perform time window segmentation based on the unified sampling rate system state dataset to obtain synchronized system state data; Step S183: Initialize the fiber optic sensor with wavelength demodulation based on the synchronization system status data to obtain the reference wavelength of the fiber optic sensor; perform spectral feature detection on the metal vapor plume based on the reference wavelength of the fiber optic sensor to obtain the plume spectral absorption data; Step S184: Perform vapor density inversion based on plume spectral absorption data to obtain vapor plume interference characteristic data; Step S185: Reconstruct the temperature field of the thermocouple array based on the synchronous system state data to obtain the servo motor temperature field distribution data; Step S186: Calculate the temperature gradient of the motor windings based on the servo motor temperature field distribution data to obtain the servo motor thermal characteristic data.

4. The method for automatic detection of control model parameters of a closed-loop CNC system in laser cutting according to claim 1, characterized in that, Step S2 involves extracting interference features from the steam plume interference data, including: Frequency domain transformation preprocessing is performed based on steam plume interference characteristic data to obtain plume frequency domain transformation preprocessing data. The frequency domain transformation preprocessing includes windowing the steam plume interference characteristic data using a Hamming window function, with the window length set to 1024 sampling points and an overlap rate of 75%. The data length is extended to 2048 points using zero-filling technology. Fast Fourier Transform is performed on the preprocessed data of the plume frequency domain transform to obtain the frequency domain data of the position feedback signal; Amplitude and phase spectra are separated based on the frequency domain data of the position feedback signal to obtain separated frequency domain feature data; Power spectral density is estimated based on the separated frequency domain characteristic data to obtain the power spectral data of the position feedback signal; Based on the power spectrum data of the location feedback signal, spectral peak detection and location are performed to obtain frequency domain characteristic data of plume interference.

5. The method for automatic detection of control model parameters of a closed-loop CNC system in laser cutting according to claim 1, characterized in that, Step S2 involves performing motor parameter drift analysis on the servo motor thermal characteristic data, including: Wavelet basis function selection is performed based on servo motor thermal characteristic data to obtain wavelet transform parameters. The wavelet basis function selection includes selecting Daubechies4 wavelet as mother wavelet function, setting the number of decomposition layers to 6, the scale parameter a varying from 1 to 64, and the step size of the translation parameter b set to 1 / 10 of the sampling interval. Continuous wavelet transform is performed on the thermal characteristic data of the servo motor according to the wavelet transform parameters to obtain the time-frequency domain data of the motor temperature. Electromagnetic constant attenuation is identified based on motor temperature time-frequency domain data to obtain raw data of electromagnetic constant attenuation. Motor thermal characteristics are identified based on the original data of electromagnetic constant decay to obtain instantaneous parameters of motor thermal characteristics; Huang transform decomposition is performed on the instantaneous parameters of motor thermal characteristics to obtain the intrinsic mode function data of thermal characteristics. Instantaneous frequency calculations are performed based on thermal characteristic intrinsic mode function data to obtain motor parameter drift characteristic data.

6. The method for automatic detection of control model parameters of a closed-loop CNC system in laser cutting according to claim 1, characterized in that, Step S2, which involves constructing a system parameter drift indicator library based on plume interference frequency domain indicator data and motor parameter drift indicator data, includes: Construct a plume interference feature vector based on plume interference frequency domain characteristic data; Construct a motor drift feature vector based on motor parameter drift symptom data; Multidimensional feature space mapping is performed based on the plume interference feature vector and the motor drift feature vector to obtain the fused system state feature space data; Clustering and pattern recognition are performed on the fusion system state feature space data to obtain a system parameter drift sign library.

7. The method for automatic detection of control model parameters of a closed-loop CNC system in laser cutting according to claim 1, characterized in that, Step S3 includes the following steps: Step S31: Initialize the system parameters using the recursive least squares method based on the system parameter drift symptom library to obtain the initial parameters of the recursive algorithm; Step S32: Real-time data acquisition is performed on the closed-loop CNC system in laser cutting to obtain system input and output data; regression matrix is ​​constructed on the system input and output data based on the initial parameters of the recursive algorithm to obtain the regression matrix of the closed-loop CNC system; Step S33: Perform recursive least squares parameter updates based on the regression matrix of the closed-loop CNC system to obtain real-time CNC system parameter estimation data; Step S34: Extract the feedforward compensation coefficients based on the parameter estimation data of the real-time CNC system to obtain the original data of the feedforward compensation coefficients; Step S35: Calculate the forgetting factor based on the original data of the feedforward compensation coefficient to obtain the forgetting factor data of the CNC system; Step S36: Adjust the weights of the original data of the feedforward compensation coefficients based on the forgetting factor data of the CNC system to obtain the weighted feedforward compensation data of the CNC system; Step S37: Perform parameter oscillation suppression on the feedforward compensation data of the weighted CNC system to obtain stable feedforward compensation parameters.

8. The method for automatic detection of control model parameters of a closed-loop CNC system in laser cutting according to claim 7, characterized in that, Step S37 includes the following steps: Step S371: Extract historical data from the sliding window based on the weighted CNC system feedforward compensation data to obtain a set of historical CNC system data windows; Step S372: Perform weighted least squares fitting based on the historical CNC system data window set to obtain the feedforward compensation coefficients of the CNC system; Step S373: Calculate the parameter change rate based on the feedforward compensation coefficient of the CNC system to obtain the compensation coefficient change rate data; Step S374: Perform threshold detection and amplitude limiting based on the rate of change of the compensation coefficient to obtain the corrected feedforward compensation coefficient; Step S375: Perform system stability assessment based on the modified feedforward compensation coefficient to obtain system stability assessment data; select an adjustment strategy based on the system stability assessment data and the preset adjustment strategy set to obtain the effective strategy adjustment parameters of the system; Step S376: Adjust the parameters based on the effective strategy of the system to perform parameter convergence verification and obtain the convergence verification data of the closed-loop CNC system; Step S377: Confirm the final parameters based on the convergence verification data of the closed-loop CNC system to obtain stable feedforward compensation parameters.

9. The method for automatic detection of control model parameters of a closed-loop CNC system in laser cutting according to claim 1, characterized in that, Step S4 includes the following steps: Step S41: Construct a mathematical model of the closed-loop system based on the stable feedforward compensation parameters to obtain the system state-space model; Step S42: Extract the characteristic polynomial coefficients based on the system state-space model to obtain the characteristic polynomial parameters of the closed-loop system; Step S43: Identify the system pole locations based on the characteristic polynomial parameters of the closed-loop system to obtain the closed-loop system pole data; Step S44: Identify stability criteria based on the closed-loop system pole data to obtain the stability coefficients of the closed-loop system; Step S45: Evaluate the stability margin of the closed-loop CNC system in laser cutting based on the stability coefficient of the closed-loop system to obtain the stability margin of the closed-loop CNC system. Step S46: Iteratively optimize the speed loop gain parameters based on the stability margin of the closed-loop CNC system to obtain the optimized speed loop control parameters.

10. The method for automatic detection of control model parameters of a closed-loop CNC system in laser cutting according to claim 1, characterized in that, Step S5 includes the following steps: Step S51: Acquire the original trajectory data of the closed-loop CNC system in laser cutting to obtain the original cutting trajectory data; Based on the speed loop optimization control parameters, reconstruct the original cutting trajectory data to obtain the standard cutting trajectory data. Step S52: Calculate the curvature based on the standard cutting trajectory data to obtain the cutting trajectory curvature distribution data; Step S53: Extract geometric feature parameters based on the curvature distribution data of the cutting trajectory to obtain the geometric feature data of the cutting trajectory; Step S54: Calculate the rate of change of tangent direction based on the geometric feature data of the cutting trajectory to obtain the feature data of the change of cutting direction; Step S55: Perform velocity planning and prediction based on the cutting direction change feature data to obtain the cutting predicted velocity distribution data; Step S56: Real-time trajectory data acquisition is performed on the closed-loop CNC system during laser cutting to obtain the current cutting trajectory data; Step S57: Perform predictive compensation control on the current cutting trajectory data based on the cutting prediction speed distribution data to obtain the final cutting trajectory control data.

Images