High-precision Z-axis module integrated driving control method and system based on force feedback

Abstract

The invention provides a high-precision Z-axis module integrated driving control method and system based on force feedback, and relates to the technical field of precise motion control, and the method comprises the steps: obtaining a position and force signal to construct a multi-dimensional state sequence, calculating a coupling relation to carry out channel decoupling, predicting a position and force evolution trajectory, constructing a bidirectional causal chain to carry out constraint adjustment, and carrying out precise motion control. According to the invention, a comprehensive driving instruction is obtained and stability monitoring adjustment is carried out, so that high-precision position and force control of the Z-axis module is realized, the response speed and stability are improved, and the adaptive capacity under complex working conditions is enhanced.

Description

TECHNICAL FIELD

[0001] The present application relates to the technical field of precision motion control, and in particular to a high-precision Z-axis module integrated driving control method and system based on force feedback. BACKGROUND

[0002] In the fields of modern precision manufacturing, semiconductors, medical devices, optical processing, etc., the high-precision control capability of the Z-axis module directly affects the product quality and production efficiency. The traditional Z-axis module control system usually adopts single position feedback control or simple force control method, which cannot meet the application scenarios with high precision, high response speed and stability requirements. With the development of industrial automation and intelligent manufacturing, higher requirements are put forward for the precision, stability and dynamic response of the Z-axis module. However, in the existing Z-axis module control technology, position control and force control are usually designed as independent systems, lacking effective coupling processing mechanism, leading to mutual interference between the two control systems under high precision requirements, making it difficult to realize accurate position-force collaborative control, lacking adaptability to system dynamic characteristic changes and external disturbances, and being difficult to maintain stable control precision under complex working conditions, lacking prediction and feedforward compensation ability of system state, relying only on feedback error correction mechanism leading to response lag, being difficult to handle dynamic mechanical effects under high-speed motion, limiting the performance of the Z-axis module in high dynamic application scenarios. SUMMARY

[0003] The embodiments of the present application provide a high-precision Z-axis module integrated driving control method and system based on force feedback, which can at least solve some problems existing in the prior art.

[0004] In a first aspect of the embodiments of the present application, a high-precision Z-axis module integrated driving control method based on force feedback is provided, comprising: obtaining position signals and force sensing signals of the Z-axis module in the motion process and synchronously aligning them according to time sequence, and calculating a multi-dimensional state sequence based on the position signals and the force sensing signals; calculating the correlation degree between the dimensions of the multi-dimensional state sequence to obtain a coupling relationship and determine a coupling influence coefficient, performing linear transformation on the multi-dimensional state sequence based on the coupling influence coefficient to obtain a channel decoupling state sequence, performing sliding window segmentation on the channel decoupling state sequence and calculating local change trend, predicting position evolution trajectory and force evolution trajectory based on the local change trend and performing time sequence alignment to obtain a predicted state trajectory pair; cross-correlation analysis is performed on the pair of predicted state trajectories to identify a bidirectional coupling influence relationship between the position and the force and to construct a bidirectional causal chain, dynamic constraint adjustment is performed on the pair of predicted state trajectories based on the bidirectional causal chain to obtain a pair of corrected predicted trajectories, a target tracking deviation of the pair of corrected predicted trajectories is calculated and a channel control component is solved, and a comprehensive driving instruction is obtained by superimposing the channel control component; According to the comprehensive driving instruction, a driving response stability index of the Z-axis module is monitored, and if the stability index exceeds a preset stability threshold, amplitude limitation and rate constraint are performed on the comprehensive driving instruction to obtain a safe driving instruction and execute the safe driving instruction.

[0005] In an optional embodiment, A position signal and a force sensing signal of the Z-axis module during movement are acquired and are synchronized and aligned in time sequence, and a multi-dimensional state sequence is calculated based on the position signal and the force sensing signal, including: A position signal of the Z-axis module is collected by a position sensor, a force sensing signal of the Z-axis module is collected by a force sensor, timestamp information of the position signal and the force sensing signal is extracted respectively, and a time offset between the position signal and the force sensing signal is identified, time sequence interpolation alignment is performed based on the time offset to obtain an aligned position signal and an aligned force signal; A first-order differential operation is performed on the aligned position signal to obtain an instantaneous velocity signal, a second-order differential operation is performed on the aligned position signal to obtain an instantaneous acceleration signal, the instantaneous velocity signal and the instantaneous acceleration signal are combined as a motion derivative component, and high-frequency noise suppression processing is performed on the aligned force signal to extract an effective frequency band component as a net force component; The aligned position signal, the motion derivative component, and the net force component are vector spliced in the time dimension to obtain the multi-dimensional state sequence.

[0006] In an optional embodiment, A coupling relationship is obtained by calculating a correlation degree between each dimension component in the multi-dimensional state sequence to determine a coupling influence coefficient, and a channel decoupling state sequence is obtained by linear transformation of the multi-dimensional state sequence based on the coupling influence coefficient, including: Cross-correlation operation is performed on any two dimension components in the multi-dimensional state sequence to obtain a correlation coefficient matrix, eigenvalue decomposition is performed on the correlation coefficient matrix, a feature vector with an eigenvalue greater than a preset coupling judgment threshold is extracted as a target feature vector, and a coupled dimension combination is obtained based on the target feature vector to identify a dimension combination with a coupling relationship; A joint state space is constructed based on the position dimension component and the force dimension component in the coupled dimension combination. The joint state space is reconstructed into a coupled phase space trajectory. Dynamic time warping analysis is performed on the coupled phase space trajectory to extract the time alignment mapping relationship between the evolution path of the position dimension component and the evolution path of the force dimension component. Based on the time alignment mapping relationship, a bidirectional regression fitting is performed on the position dimension component and the force dimension component to obtain a bidirectional transfer function pair, wherein the bidirectional transfer function pair constitutes a coupling relationship. Frequency domain analysis is performed on the bidirectional transfer function pair to obtain frequency response characteristics and the gain coefficient of the dominant frequency component is extracted as the coupling influence coefficient. Based on the coupling influence coefficient and the target feature vector, a decoupling transformation matrix is ​​constructed, and a matrix left multiplication transformation is performed on the multidimensional state sequence to obtain the position channel state sequence and the force channel state sequence, which are then combined to obtain the channel decoupling state sequence.

[0007] In one alternative implementation, The channel decoupling state sequence is segmented by a sliding window and local change trends are calculated. Based on these local change trends, position evolution trajectories and force evolution trajectories are predicted and time-aligned to obtain predicted state trajectory pairs, including: A sliding window is set for the channel decoupling state sequence. The channel decoupling state sequence is divided into multiple window segments according to a preset step size. Multidimensional linear fitting is performed on the data in each window segment. The fitting slope of the position channel dimension is calculated to obtain the position change rate, and the fitting slope of the force channel dimension is calculated to obtain the force change rate. The position change rate and the force change rate are combined to obtain the local change trend. Based on the local change trend, the end value of the position channel in the channel decoupling state sequence is extracted as the position starting point. The position evolution trajectory is obtained by linear extrapolation from the position starting point with the position change rate as the extension rate. The end value of the force channel in the channel decoupling state sequence is extracted as the force starting point. The force evolution trajectory is obtained by linear extrapolation from the force starting point with the force change rate as the extension rate. Extract the first timestamp sequence corresponding to the position evolution trajectory and the second timestamp sequence corresponding to the force evolution trajectory. Perform time-series alignment on the first timestamp sequence and the second timestamp sequence to obtain a unified time reference. Based on the unified time reference, map the position evolution trajectory and the force evolution trajectory to the same time coordinate system to obtain a predicted state trajectory pair.

[0008] In one alternative implementation, Cross-correlation analysis is performed on the predicted state trajectory pairs to identify the bidirectional coupling influence relationship between position and force and to construct a bidirectional causal chain. Based on the bidirectional causal chain, the predicted state trajectory pairs are dynamically constrained and adjusted to obtain corrected predicted trajectory pairs, including: From the predicted state trajectory pair, extract the position numerical sequence corresponding to the position evolution trajectory and the force numerical sequence corresponding to the force evolution trajectory. Differentiate the position numerical sequence to obtain the velocity sequence. Differentiate the velocity sequence to obtain the acceleration sequence. Calculate the cross-correlation function between the acceleration sequence and the force numerical sequence. At the peak time of the cross-correlation function, extract the reciprocal of the dynamic response delay of the force to the position to obtain the force causality intensity. Integrate the force numerical sequence to obtain the impulse sequence. Calculate the partial cross-information between the impulse sequence and the position numerical sequence and extract the energy transfer efficiency of the position to the force to obtain the position causality intensity. Based on the force causality intensity and the position causality intensity, identify the bidirectional coupling influence relationship. Position and force states are set as nodes in a causal graph. Directed connections from position nodes to force nodes and from force nodes to position nodes are established based on bidirectional coupling influence relationships. The bidirectional causal chain is obtained by combining the directed connections. Path analysis is performed on the bidirectional causal chain to identify the forward propagation path from the position node to the force node and the backward propagation path from the force node to the position node, and a bidirectional coupling constraint relationship is established. Based on the bidirectional coupling constraint relationship, the predicted state trajectory pair is corrected to obtain the corrected predicted trajectory pair.

[0009] In one alternative implementation, Calculate the target tracking deviation of the corrected predicted trajectory pair and solve for the channel control components. Superimpose the channel control components to obtain the comprehensive drive command, including: The current task corresponding to the corrected predicted trajectory pair is identified by feature analysis. The target trajectory corresponding to the current task is extracted from the preset trajectory library. The corrected predicted trajectory pair and the target trajectory are embedded into a low-dimensional manifold space by manifold learning algorithm and geodesic distance is calculated in the low-dimensional manifold space. The dominant direction of trajectory deviation is extracted based on the geodesic distance. The target tracking deviation is determined based on the geodesic distance and the dominant direction. The membership value corresponding to the target tracking deviation is calculated based on the preset Gaussian membership function. The target tracking deviation is mapped to a fuzzy linguistic variable based on the membership value. Fuzzy inference is performed on the fuzzy linguistic variable based on the preset fuzzy rule base to obtain the activation intensity of each fuzzy rule. Fuzzy control output is calculated based on the activation intensity and the fuzzy rule. The fuzzy control output is defuzzified using the centroid method to obtain the channel control component. The comprehensive driving command is obtained by superimposing the amplitude normalized components of the channel control components.

[0010] In one alternative implementation, The Z-axis module is driven according to the integrated drive command, and the stability index of the drive response is monitored. If the stability index exceeds a preset stability threshold, the amplitude and rate of change of the integrated drive command are limited to obtain a safe drive command and executed, including: The integrated drive command is converted into a pulse width modulation signal and applied to the drive motor of the Z-axis module. Real-time force response signal and real-time position response signal are collected. The real-time force response signal and real-time position response signal are decomposed into multiple layers by wavelet packet decomposition algorithm to extract energy value and construct energy feature vector. The center of the minimum enclosing hypersphere of the energy feature vector is calculated by support vector data description algorithm. The Mahalanobis distance from the energy feature vector to the center of the sphere is calculated to obtain the stability index. When the stability index is greater than the preset stability threshold, the amplitude of the integrated drive command is extracted and compared with the preset amplitude upper limit. The amplitude exceeding the preset amplitude upper limit is clipped to the preset amplitude upper limit to obtain the amplitude clipped command. The difference between the amplitude clipped command at adjacent time points is calculated to obtain the amplitude change rate. The amplitude change rate is compared with the preset change rate upper limit and the part exceeding the preset change rate upper limit is smoothed to obtain the safety drive command. The safety drive command is converted into a pulse width modulation signal and applied to the drive motor.

[0011] A second aspect of the present invention provides a high-precision Z-axis module integrated drive control system based on force feedback, comprising: The first unit is used to acquire the position signal and force sensing signal of the Z-axis module during the motion process and synchronize them according to the time sequence, and calculate a multi-dimensional state sequence based on the position signal and the force sensing signal. The second unit is used to calculate the degree of correlation between the components of each dimension in the multidimensional state sequence to obtain the coupling relationship and determine the coupling influence coefficient. Based on the coupling influence coefficient, the multidimensional state sequence is linearly transformed to obtain the channel decoupled state sequence. The channel decoupled state sequence is segmented by a sliding window and the local change trend is calculated. Based on the local change trend, the position evolution trajectory and the force evolution trajectory are predicted and time-series aligned to obtain the predicted state trajectory pair. The third unit is used to perform cross-correlation analysis on the predicted state trajectory pair, identify the bidirectional coupling influence relationship between position and force and construct a bidirectional causal chain, dynamically adjust the predicted state trajectory pair based on the bidirectional causal chain to obtain a corrected predicted trajectory pair, calculate the target tracking deviation of the corrected predicted trajectory pair and solve for the channel control component, and superimpose the channel control component to obtain a comprehensive driving command. The fourth unit is used to apply drive to the Z-axis module according to the comprehensive drive command and monitor the stability index of the drive response. If the stability index exceeds the preset stability threshold, the comprehensive drive command is subjected to amplitude limitation and rate of change constraint to obtain a safe drive command and execute it.

[0012] A third aspect of the present invention provides an electronic device, comprising: A processor and a memory for storing processor-executable instructions, wherein the processor is configured to invoke instructions stored in the memory to perform the aforementioned method.

[0013] A fourth aspect of the present invention provides a computer-readable storage medium having stored thereon computer program instructions that, when executed by a processor, implement the aforementioned method.

[0014] In this invention, a high-precision Z-axis module integrated drive control method based on force feedback is used to achieve high-precision synchronous processing and channel decoupling of position and force signals, effectively eliminating interference between signals, improving the accuracy of state estimation and the sensitivity of system response. By constructing a bidirectional causal chain and dynamic constraint adjustment, accurate modeling and control of the complex interaction between position and force during the movement of the Z-axis module are achieved, effectively solving the limitations of traditional control methods in dealing with coupled systems, significantly improving the system's adaptability under complex working conditions, and introducing a stability monitoring and safety drive mechanism, which can evaluate the system state in real time and make necessary constraint adjustments to ensure the system stability and safety during high-precision drive. Attached Figure Description

[0015] Figure 1 This is a flowchart illustrating the high-precision Z-axis module integrated drive control method based on force feedback according to an embodiment of the present invention. Figure 2 This is a flowchart illustrating the bidirectional causal coupling analysis and trajectory correction of the high-precision Z-axis module integrated drive control method based on force feedback, as described in an embodiment of the present invention. Detailed Implementation

[0016] To make the objectives, technical solutions, and advantages of the embodiments of the present invention clearer, the technical solutions of the embodiments of the present invention will be clearly and completely described below with reference to the accompanying drawings. Obviously, the described embodiments are only some embodiments of the present invention, and not all embodiments. Based on the embodiments of the present invention, all other embodiments obtained by those skilled in the art without creative effort are within the scope of protection of the present invention.

[0017] The technical solution of the present invention will be described in detail below with reference to specific embodiments. These specific embodiments can be combined with each other, and the same or similar concepts or processes may not be described again in some embodiments.

[0018] Figure 1 This is a flowchart illustrating the high-precision Z-axis module integrated drive control method based on force feedback according to an embodiment of the present invention. Figure 1 As shown, the method includes: The position signal and force sensing signal of the Z-axis module during the motion process are acquired and synchronized according to the time sequence. A multi-dimensional state sequence is calculated based on the position signal and the force sensing signal. The degree of correlation between the components of each dimension in the multidimensional state sequence is calculated to obtain the coupling relationship and the coupling influence coefficient is determined. Based on the coupling influence coefficient, the multidimensional state sequence is linearly transformed to obtain the channel decoupled state sequence. The channel decoupled state sequence is segmented by a sliding window and the local change trend is calculated. Based on the local change trend, the position evolution trajectory and the force evolution trajectory are predicted and time-series aligned to obtain the predicted state trajectory pair. Cross-correlation analysis is performed on the predicted state trajectory pair to identify the bidirectional coupling influence relationship between position and force and construct a bidirectional causal chain. Based on the bidirectional causal chain, dynamic constraint adjustment is performed on the predicted state trajectory pair to obtain a corrected predicted trajectory pair. The target tracking deviation of the corrected predicted trajectory pair is calculated and the channel control component is obtained. The channel control component is superimposed to obtain the comprehensive drive command. The Z-axis module is driven according to the integrated drive command, and the stability index of the drive response is monitored. If the stability index exceeds the preset stability threshold, the amplitude limit and rate of change constraint of the integrated drive command are applied to obtain a safe drive command and execute it.

[0019] In one alternative implementation, The position signal and force sensing signal of the Z-axis module during the motion process are acquired and synchronized according to the time sequence. Based on the position signal and the force sensing signal, a multidimensional state sequence is calculated, including: The Z-axis module's position signal is acquired by a position sensor, and the Z-axis module's force sensing signal is acquired by a force sensor. The timestamp information of the position signal and the force sensing signal is extracted respectively, and the time offset between the position signal and the force sensing signal is identified. Based on the time offset, timing interpolation and alignment are performed to obtain the aligned position signal and the aligned force signal. The instantaneous velocity signal is obtained by performing a first-order differential operation on the alignment position signal, and the instantaneous acceleration signal is obtained by performing a second-order differential operation on the alignment position signal. The instantaneous velocity signal and the instantaneous acceleration signal are combined into a motion derivative component. The alignment force signal is subjected to high-frequency noise suppression processing, and the effective frequency band component is extracted and used as the net force component. The alignment position signal, the motion derivative component, and the net force component are concatenated as vectors along the time dimension to obtain the multidimensional state sequence.

[0020] The Z-axis module's position signal is acquired via a position sensor, and its force signal is acquired via a force sensor. The position sensor can be a photoelectric encoder, Hall effect sensor, or electromagnetic displacement sensor, with a sampling frequency set at 1000Hz and an accuracy of 0.01 mm. The force sensor can be a piezoelectric or strain gauge force sensor, also with a sampling frequency set at 1000Hz, a measurement range of 0-100 Newtons, and an accuracy of 0.05 Newtons.

[0021] After signal acquisition, timestamps are extracted from the position and force sensing signals. Each acquired data point contains a timestamp indicating the precise moment of data acquisition. Due to differences in hardware and signal transmission delays, there is usually a time offset between the position and force sensing signals. The average time offset is calculated by comparing the timestamp sequences of the two signals. In practical applications, this time offset is typically between 5 and 20 milliseconds.

[0022] Based on the calculated time offset, temporal interpolation alignment is performed. Cubic spline interpolation is used to interpolate signals with lower sampling rates to the same time points as signals with higher sampling rates. For example, if the position signal timestamps are 0, 10, 20, and 30 milliseconds, the force signal timestamps are 3, 13, 23, and 33 milliseconds, and the time offset is 3 milliseconds, then the force signal is aligned to the time points of 0, 10, 20, and 30 milliseconds through interpolation, resulting in aligned position and force signals at the same time points.

[0023] The obtained alignment position signal is subjected to first-order differentiation to calculate the instantaneous velocity at each time point. The calculation involves dividing the difference in position values ​​between adjacent time points by the time interval. To reduce the influence of noise, the center difference method is used for calculation, employing a weighted average of the position values ​​from multiple points before and after the center difference. For edge points, forward difference or backward difference methods are used for calculation.

[0024] The second-order differential operation is performed on the alignment position signal to obtain the instantaneous acceleration signal. The instantaneous acceleration is calculated based on the ratio of the velocity change at adjacent time points to the time interval. To improve calculation accuracy, the central difference method is also used, and a sliding window averaging method is employed to reduce the influence of noise. The window size is set to 5 sampling points, which can effectively suppress high-frequency noise while preserving signal characteristics.

[0025] The calculated instantaneous velocity and instantaneous acceleration signals are combined into a motion derivative component. During the merging process, the velocity and acceleration are normalized to ensure consistent dimensions. The normalization range for the velocity signal is set to -1 to +1, corresponding to the actual velocity range of -500 mm / s to +500 mm / s; the normalization range for the acceleration signal is also set to -1 to +1, corresponding to the actual acceleration range of -5000 mm / s² to +5000 mm / s².

[0026] High-frequency noise suppression is performed on the force signal. Force sensor signals are often affected by mechanical vibration, electromagnetic interference, and other factors, containing a large amount of high-frequency noise. A low-pass filter is used for processing, with a cutoff frequency set to 50 Hz, which can retain the effective information of the force signal while filtering out most of the high-frequency noise. The filter is a Butterworth filter with an order of 4 and a transition band width of 10 Hz.

[0027] After noise suppression, the effective frequency band component of the force signal is extracted as the net force component. Through spectrum analysis, the effective frequency band of the Z-axis module force signal is determined to be 0-30 Hz. A bandpass filter is used to extract the signal components in this frequency band. The filter parameters are designed as follows: passband 5-25 Hz, stopband 0-2 Hz and 35-500 Hz, passband ripple not exceeding 0.5 dB, and stopband attenuation greater than 40 dB.

[0028] The alignment position signal, motion derivative components, and net force components are concatenated into vectors along the time dimension to form a multidimensional state sequence. The vector concatenation method involves combining the corresponding position, velocity, acceleration, and net force values ​​into a four-dimensional vector at each time point. For data with a sampling duration of 10 seconds and a sampling frequency of 1000 Hz, the final multidimensional state sequence is a 10,000-row, 4-column matrix.

[0029] In this embodiment, by identifying and aligning the time offsets of the position and force signals through time-series interpolation, the consistency of the two types of signals in the time dimension is ensured, avoiding data distortion caused by sampling delays or asynchrony. By performing first- and second-order differential operations on the aligned position signals, more accurate motion derivative features are obtained, which can comprehensively reflect the dynamic response process of the Z-axis module. By suppressing high-frequency noise and extracting effective frequency bands for the force signals, the purity and feature stability of the signals are significantly improved. By fusing the position signals, motion derivative components, and net force components in the time dimension, a multi-dimensional state sequence is constructed, enhancing the joint perception capability of the motion and force states of the Z-axis module, and improving the reliability and response accuracy of subsequent motion control and state assessment.

[0030] In one alternative implementation, Calculating the correlation between the components of each dimension in the multidimensional state sequence to obtain the coupling relationship and determining the coupling influence coefficient, and performing a linear transformation on the multidimensional state sequence based on the coupling influence coefficient to obtain the channel decoupled state sequence includes: Cross-correlation operation is performed on any two dimensional components in the multidimensional state sequence to obtain a correlation coefficient matrix. Eigenvalue decomposition is performed on the correlation coefficient matrix, and feature vectors with eigenvalues ​​greater than a preset coupling judgment threshold are extracted as target feature vectors. Based on the target feature vectors, dimensional combinations with coupling relationships are identified to obtain coupled dimensional combinations. A joint state space is constructed based on the position dimension component and the force dimension component in the coupled dimension combination. The joint state space is reconstructed into a coupled phase space trajectory. Dynamic time warping analysis is performed on the coupled phase space trajectory to extract the time alignment mapping relationship between the evolution path of the position dimension component and the evolution path of the force dimension component. Based on the time alignment mapping relationship, a bidirectional regression fitting is performed on the position dimension component and the force dimension component to obtain a bidirectional transfer function pair, wherein the bidirectional transfer function pair constitutes a coupling relationship. Frequency domain analysis is performed on the bidirectional transfer function pair to obtain frequency response characteristics and the gain coefficient of the dominant frequency component is extracted as the coupling influence coefficient. Based on the coupling influence coefficient and the target feature vector, a decoupling transformation matrix is ​​constructed, and a matrix left multiplication transformation is performed on the multidimensional state sequence to obtain the position channel state sequence and the force channel state sequence, which are then combined to obtain the channel decoupling state sequence.

[0031] Cross-correlation is performed on any two dimensional components in a multidimensional state sequence to calculate the correlation coefficient. In a four-dimensional state sequence, which includes position, velocity, acceleration, and force components, the correlation coefficient matrix is ​​obtained by calculating the cross-correlation coefficients between any two dimensions. The cross-correlation coefficients are calculated based on the average of the normalized products of the two dimensional components. For example, the data points for the position and force components are normalized separately, the product of the normalized values ​​at the corresponding time points is calculated, and the average is taken to obtain the cross-correlation coefficient. Cross-correlation coefficients are calculated pairwise for all dimensional components, ultimately resulting in a four-by-four correlation coefficient matrix.

[0032] Eigenvalue decomposition is performed on the correlation coefficient matrix to obtain the corresponding eigenvalues ​​and eigenvectors. Eigenvalue decomposition is implemented using an iterative algorithm; for a 4x4 correlation coefficient matrix, four eigenvalues ​​and four corresponding eigenvectors are obtained. The eigenvalues ​​represent the correlation strength between different dimensions, and the eigenvectors represent the correlation pattern. A preset coupling threshold of 0.6 is set, and eigenvectors with eigenvalues ​​greater than this threshold are extracted as target eigenvectors. In practical applications, typically 1-2 target eigenvectors are obtained.

[0033] This method identifies coupled dimensional combinations based on target feature vectors. It examines the weight coefficients of each dimension in each target feature vector. If the absolute values ​​of the weight coefficients of two dimensions in the feature vector are both greater than 0.4, then these two dimensions are considered coupled, and all coupled dimensional combinations are identified. In the Z-axis module, position and force components typically have a strong coupling relationship and are therefore identified as coupled dimensional combinations.

[0034] A joint state space is constructed with position components as the x-axis and force components as the y-axis. For the sampled multidimensional state sequence, the position and force values ​​at each time point constitute a point in the joint state space. Data from all time points are plotted in the joint state space to form a trajectory. To achieve phase space reconstruction, a time delay parameter is introduced, set to 5 sampling points. The original trajectory is time-delayed, and a coupled phase space trajectory is reconstructed in the multidimensional space.

[0035] Dynamic time warping analysis is performed on the coupled phase space trajectory to extract the time alignment mapping relationship between the position component and the force component. The dynamic time warping algorithm finds the time correspondence by calculating the optimal alignment path between two time series. With a window size of 10 sampling points and a step size of 1 sampling point, the similarity between the position component and the force component at different time points is calculated, a similarity matrix is ​​constructed, and a dynamic programming algorithm is used to determine the path with the highest similarity; this path represents the optimal time alignment mapping relationship between the position component and the force component.

[0036] Based on the time-aligned mapping relationship, a bidirectional regression fitting is performed on the position and force dimension components to obtain bidirectional transfer function pairs. The transfer function from position to force describes how changes in position affect changes in position, while the transfer function from force to position describes how changes in force affect changes in position. A polynomial regression model is used for fitting, with the polynomial order set to 3. During the fitting process, the least squares method is used to determine the polynomial coefficients. The fitted bidirectional transfer functions form a coupling relationship, representing the mapping relationship from position to force and from force to position, respectively.

[0037] Frequency domain analysis is performed on the bidirectional transfer function pair to obtain its frequency response characteristics. The time-domain transfer function is converted into a frequency-domain frequency response function using a discrete Fourier transform. The amplitude and phase characteristics of the frequency response function are calculated to analyze the transmission relationship between position and force at different frequencies. The gain coefficient of the dominant frequency component is extracted as the coupling influence coefficient. The dominant frequency is typically in the 0-10 Hz range, corresponding to the main dynamic characteristics of the Z-axis module. Within the dominant frequency range, the point with the largest frequency response amplitude is identified; this amplitude is the coupling influence coefficient.

[0038] A decoupling transformation matrix is ​​constructed based on the coupling influence coefficient and the target eigenvector. The design goal of the decoupling transformation matrix is ​​to reduce the coupling relationship between the position component and the force component. The weight coefficients of each dimension of the target eigenvector are multiplied by the coupling influence coefficient to form the elements of the decoupling transformation matrix. The dimension of the decoupling transformation matrix is ​​the same as the dimension of the multidimensional state sequence, which is four by four.

[0039] A matrix left-multiplication transformation is performed on the multidimensional state sequence, that is, the decoupling transformation matrix is ​​multiplied by the state vector at each time point in the multidimensional state sequence to obtain the transformed state vector. In the transformed state vector, the position-related components (position, velocity, acceleration) form the position channel state sequence, and the force-related components form the force channel state sequence. Combining the position channel state sequence and the force channel state sequence yields the channel decoupled state sequence.

[0040] In this embodiment, by performing cross-correlation operations on the multidimensional state sequence and extracting feature vectors with feature values ​​greater than a preset threshold, strongly correlated signal dimensions can be identified in the high-dimensional feature space, improving the accuracy and noise resistance of coupling detection. By constructing a joint state space based on the combination of coupling dimensions and reconstructing the phase space, the nonlinear dynamic relationship between position and force is fully preserved, avoiding the neglect of complex coupling features by traditional linear analysis methods. By extracting the time mapping relationship between the two types of signals through dynamic time warping analysis, dynamic alignment is achieved under the condition of time delay and rate difference, improving the accuracy and stability of timing matching. By constructing a decoupling transformation matrix based on the coupling influence coefficient and performing channel separation on the multidimensional state sequence, the cross-interference between channels is effectively reduced, significantly improving the stability and reliability of the system in dynamic control and state analysis.

[0041] In one alternative implementation, The channel decoupling state sequence is segmented by a sliding window and local change trends are calculated. Based on these local change trends, position evolution trajectories and force evolution trajectories are predicted and time-aligned to obtain predicted state trajectory pairs, including: A sliding window is set for the channel decoupling state sequence. The channel decoupling state sequence is divided into multiple window segments according to a preset step size. Multidimensional linear fitting is performed on the data in each window segment. The fitting slope of the position channel dimension is calculated to obtain the position change rate, and the fitting slope of the force channel dimension is calculated to obtain the force change rate. The position change rate and the force change rate are combined to obtain the local change trend. Based on the local change trend, the end value of the position channel in the channel decoupling state sequence is extracted as the position starting point. The position evolution trajectory is obtained by linear extrapolation from the position starting point with the position change rate as the extension rate. The end value of the force channel in the channel decoupling state sequence is extracted as the force starting point. The force evolution trajectory is obtained by linear extrapolation from the force starting point with the force change rate as the extension rate. Extract the first timestamp sequence corresponding to the position evolution trajectory and the second timestamp sequence corresponding to the force evolution trajectory. Perform time-series alignment on the first timestamp sequence and the second timestamp sequence to obtain a unified time reference. Based on the unified time reference, map the position evolution trajectory and the force evolution trajectory to the same time coordinate system to obtain a predicted state trajectory pair.

[0042] A sliding window is set for the channel decoupling state sequence to capture short-term data change trends. The window size is set to 200 sampling points, equivalent to a 200-millisecond data segment at a sampling frequency of 1000Hz. The channel decoupling state sequence is segmented by sliding according to a preset step size, resulting in multiple window segments. The step size is set to 20 sampling points, i.e., each slide is 20 milliseconds. For a 10-second channel decoupling state sequence, approximately 490 window segments can be obtained.

[0043] Multidimensional linear fitting is performed on the data within each window segment. For the position channel dimension, which includes three components: position, velocity, and acceleration, the least squares method is used for linear fitting. During the fitting process, the time point within the window is used as the independent variable, and the values ​​of each component are used as the dependent variable to obtain the fitted line. The fitting slope of the position channel dimension is calculated to obtain the rate of change of position. The rate of change of position contains three values, corresponding to the change trends of position, velocity, and acceleration, respectively. Linear fitting is performed on the force channel dimension, and the fitting slope of the force channel dimension is calculated to obtain the rate of change of force. The force channel includes one component: force value, and the rate of change of force is a scalar value. Combining the rate of change of position and the rate of change of force yields the local change trend. The local change trend is a four-dimensional vector containing the rates of change of four components: position, velocity, acceleration, and force.

[0044] Based on local change trends, the end values ​​of the position channels in the channel decoupling state sequence are extracted as the position starting point. The position starting point is the position, velocity, and acceleration values ​​corresponding to the last time point of the window. Using the position change rate as the extension rate, linear extrapolation is performed starting from the position starting point to obtain the position evolution trajectory. The linear extrapolation time length is set to 50 milliseconds, equivalent to 50 sampling points. During the extrapolation process, for each prediction time point, the predicted value for the current time point is obtained by adding the position change rate multiplied by the time increment to the value of the position starting point. For example, if the position value of the position starting point is 50 millimeters and the position change rate is 2 millimeters per second, then the predicted position value after 10 milliseconds is 50.02 millimeters. The same method is used to extrapolate the velocity and acceleration values.

[0045] The final value of the force channel in the decoupled state sequence is extracted as the force initiation point. The force initiation point is the force value corresponding to the last time point of the window. Using the force change rate as the extension rate, linear extrapolation is performed starting from the force initiation point to obtain the force evolution trajectory. The extrapolation time is set to 50 milliseconds. During the extrapolation process, for each prediction time point, the predicted force value at that time point is obtained by adding the force change rate and the time increment to the force initiation point value. For example, if the force value at the force initiation point is 10 Newtons and the force change rate is 0.5 Newtons per second, then the predicted force value after 10 milliseconds is 10.005 Newtons.

[0046] Extract the first timestamp sequence corresponding to the position evolution trajectory. The first timestamp sequence starts from the last time point of the window and increments by 1 millisecond sampling intervals, containing a total of 50 time points. Similarly, extract the second timestamp sequence corresponding to the force evolution trajectory. Since the position and force evolution trajectories use the same prediction duration and sampling interval, the two timestamp sequences are usually consistent. However, in some cases, slight deviations may exist due to differences in computational precision or hardware.

[0047] The first and second timestamp sequences are time-aligned to obtain a unified time reference. The time alignment uses nearest-neighbor interpolation to ensure that each point in the two timestamp sequences corresponds to a point on the unified time reference. The time interval of the unified time reference is set to 1 millisecond, consistent with the original sampling interval. Based on the unified time reference, the position evolution trajectory and force evolution trajectory are mapped to the same time coordinate system to obtain predicted state trajectory pairs. Each predicted state trajectory pair contains four components: position, velocity, acceleration, and force, with each component having a corresponding predicted value at the same time point.

[0048] When the Z-axis module performs precision machining tasks, the aforementioned method can predict the position and force change trends in the near future. For example, when performing a pressure detection task, the Z-axis module moves downwards at a speed of 5 mm / s. Upon contact with the workpiece surface, the force value rapidly increases from 0 Newtons. Using this method, the force change trend can be detected 5-10 milliseconds in advance, allowing for corresponding adjustments to the position control strategy to avoid excessive impact force.

[0049] In this embodiment, by setting a sliding window and performing multidimensional linear fitting on the channel decoupled state sequence, local data features can be captured in a piecewise form. This suppresses instantaneous noise interference in the rate of change calculation, improving the smoothness and stability of trend estimation. By calculating the rate of change of position and the rate of change of force separately within the window and combining them into a local trend, the local dynamic coupling direction of position and force can be reflected simultaneously, significantly improving the ability to distinguish time-varying features. Linear extrapolation based on the local trend allows for rapid deduction of future states without relying on complex prediction models, improving the real-time performance and response speed of trajectory prediction. In one alternative implementation, Cross-correlation analysis is performed on the predicted state trajectory pairs to identify the bidirectional coupling influence relationship between position and force and to construct a bidirectional causal chain. Based on the bidirectional causal chain, the predicted state trajectory pairs are dynamically constrained and adjusted to obtain corrected predicted trajectory pairs, including: From the predicted state trajectory pair, extract the position numerical sequence corresponding to the position evolution trajectory and the force numerical sequence corresponding to the force evolution trajectory. Differentiate the position numerical sequence to obtain the velocity sequence. Differentiate the velocity sequence to obtain the acceleration sequence. Calculate the cross-correlation function between the acceleration sequence and the force numerical sequence. At the peak time of the cross-correlation function, extract the reciprocal of the dynamic response delay of the force to the position to obtain the force causality intensity. Integrate the force numerical sequence to obtain the impulse sequence. Calculate the partial cross-information between the impulse sequence and the position numerical sequence and extract the energy transfer efficiency of the position to the force to obtain the position causality intensity. Based on the force causality intensity and the position causality intensity, identify the bidirectional coupling influence relationship. Position and force states are set as nodes in a causal graph. Directed connections from position nodes to force nodes and from force nodes to position nodes are established based on bidirectional coupling influence relationships. The bidirectional causal chain is obtained by combining the directed connections. Path analysis is performed on the bidirectional causal chain to identify the forward propagation path from the position node to the force node and the backward propagation path from the force node to the position node, and a bidirectional coupling constraint relationship is established. Based on the bidirectional coupling constraint relationship, the predicted state trajectory pair is corrected to obtain the corrected predicted trajectory pair.

[0050] The position value sequence corresponding to the position evolution trajectory and the force value sequence corresponding to the force evolution trajectory are extracted from the predicted state trajectory pairs. The position value sequence is a vector of length 50, with each element corresponding to a position value in the predicted trajectory. The force value sequence is also a vector of length 50, with each element corresponding to a force value in the predicted trajectory. The velocity sequence is obtained by differentiating the position value sequences. The differentiation operation is performed by dividing the difference between the position values ​​of two adjacent time points by the time interval. For example, if the position values ​​of two adjacent time points are 10 mm and 10.05 mm, and the time interval is 1 ms, then the velocity value at that time point is 50 mm / s. To reduce the influence of noise, the central difference method is used to calculate the differentiation, that is, using the position values ​​of two consecutive points to calculate the velocity of the current point.

[0051] The velocity sequence is differentiated to obtain the acceleration sequence. Using the same differentiation method as described above, the difference between the velocity values ​​at two adjacent time points is calculated and divided by the time interval to obtain the acceleration value. To further improve the calculation accuracy, a low-pass filter is applied to the acceleration sequence for smoothing, with the filter cutoff frequency set to 100 Hz.

[0052] Calculate the cross-correlation function between the acceleration sequence and the force numerical sequence. The cross-correlation function represents the similarity between the two sequences at different time delays. During the calculation, both the acceleration and force numerical sequences are normalized to have a mean of 0 and a standard deviation of 1. For different time delay values, the average of the element-wise products of the two sequences is calculated to obtain the cross-correlation function. The time delay range is set to -10 to +10 milliseconds, with a step size of 1 millisecond. A total of 21 cross-correlation values ​​are calculated.

[0053] The force causality strength is obtained by extracting the reciprocal of the dynamic response delay of the force to the position at the peak of the cross-correlation function. The peak of the cross-correlation function represents the time delay at which the two sequences are most similar. Find the maximum value of the cross-correlation function and its corresponding time delay. In the Z-axis module, a typical peak time delay is 3-5 milliseconds, indicating that the change in force lags behind the change in acceleration by 3-5 milliseconds. The force causality strength is calculated as the reciprocal of the peak time delay; for example, if the peak delay is 4 milliseconds, the force causality strength is 0.25.

[0054] The impulse sequence is obtained by integrating the force value sequence. The integration operation is achieved by calculating and accumulating the product of the force value and the time interval. For data with a sampling interval of 1 millisecond, the impulse increment at each time point is equal to the force value at that time point multiplied by 0.001 seconds. Each element in the impulse sequence represents the cumulative effect of the force from the start time to the current time. For example, if the force value at a certain moment is constant at 10 Newtons and lasts for 5 milliseconds, then the impulse value at that moment is 0.05 Newton-seconds.

[0055] The partial mutual information (PMI) of the impulse and position numerical sequences is calculated, and the energy transfer efficiency from position to force is extracted to obtain the positional causality strength. PMI measures the nonlinear dependency between the two sequences, considering the influence of other variables. During calculation, the probability distributions of the impulse and position sequences are estimated using kernel density estimation with a bandwidth of 0.1. The expected value of the logarithmic ratio of the product of the joint distribution and the marginal distributions is calculated to obtain the PMI value. The energy transfer efficiency from position to force is defined as the PMI value divided by the theoretical maximum MTI, typically between 0 and 1. The positional causality strength is this energy transfer efficiency, with typical values ​​between 0.4 and 0.7.

[0056] The bidirectional coupling influence relationship is identified based on the force causality strength and the position causality strength. This relationship describes the effect of position change on force and the effect of force change on position. If the force causality strength is greater than the position causality strength, the force has a stronger influence on position; conversely, the position has a stronger influence on force. When the Z-axis module performs precision pressure control tasks, the typical force causality strength is 0.3, and the position causality strength is 0.5, indicating that position change has a more significant effect on force.

[0057] Position and force states are set as nodes in a causal graph. Directed connections are established from position nodes to force nodes and from force nodes to position nodes based on bidirectional coupling influence relationships. The connection weight from a position node to a force node is set as the position causal strength, and the connection weight from a force node to a position node is set as the force causal strength. For example, the connection weight from position to force is 0.5, and the connection weight from force to position is 0.3. Combining these directed connections yields bidirectional causal chains, forming a closed-loop causal relationship network between position and force.

[0058] Path analysis is performed on the bidirectional causal chain to identify the forward propagation path from the position node to the force node and the backward propagation path from the force node to the position node. The path analysis employs a depth-first search algorithm, starting from the initial node and exploring all possible paths along directed connections until the target node is reached. In this example, the forward propagation path is position → force, and the backward propagation path is force → position, forming a complete closed-loop propagation path: position → force → position.

[0059] A two-way coupling constraint relationship is established based on the identified propagation path. This relationship describes the dynamic equilibrium condition between position changes and force changes. According to Newton's laws of mechanics, the second derivative of position (acceleration) is proportional to the applied force, with the proportionality constant being the reciprocal of the mass. In the Z-axis module, the mass is typically 2-5 kg. The two-way coupling constraint relationship is expressed as follows: acceleration should equal force divided by mass, and the rate of change of force should have a certain proportional relationship with the rate of change of position.

[0060] The predicted trajectory pairs are corrected based on the bidirectional coupling constraint relationship to obtain corrected predicted trajectory pairs. During the correction process, the theoretical acceleration value is calculated based on the predicted force value and mass, and compared with the predicted acceleration value to calculate the deviation. If the deviation exceeds a preset threshold (typically 0.1 m / s²), the predicted values ​​of position, velocity, acceleration, and force are adjusted according to the bidirectional coupling constraint relationship to satisfy the physical constraints. For example, if the predicted force value at a certain moment is 20 Newtons and the mass is 4 kg, the theoretical acceleration should be 5 m / s²; if the predicted acceleration is 4.7 m / s², the deviation is 0.3 m / s², exceeding the threshold, and correction is required.

[0061] In this embodiment, by performing differential and integral analysis on the position and force numerical sequences in the predicted state trajectory pair, interaction features can be extracted at both the time and energy domains, improving the physical consistency and interpretability of causal relationship identification. By calculating energy transfer efficiency based on the partial mutual information of impulse and position, the feedback effect of position on force is effectively quantified, extending the causal relationship judgment from unidirectional dependence to bidirectional energy interaction, thus improving the identification accuracy of complex interaction mechanisms. By establishing a bidirectional causal chain structure of position and force states, a structured expression of causal direction is realized, which facilitates the differentiation between forward propagation and reverse feedback effects in path analysis. The predicted trajectory is corrected based on the bidirectional coupling constraint relationship, so that the prediction results can take into account both physical consistency and temporal rationality, reducing the cumulative deviation under strong coupling conditions.

[0062] Figure 2 This is a flowchart illustrating the bidirectional causal coupling analysis and trajectory correction of the high-precision Z-axis module integrated drive control method based on force feedback, as described in an embodiment of the present invention.

[0063] In one alternative implementation, Calculate the target tracking deviation of the corrected predicted trajectory pair and solve for the channel control components. Superimpose the channel control components to obtain the comprehensive drive command, including: The current task corresponding to the corrected predicted trajectory pair is identified by feature analysis. The target trajectory corresponding to the current task is extracted from the preset trajectory library. The corrected predicted trajectory pair and the target trajectory are embedded into a low-dimensional manifold space by manifold learning algorithm and geodesic distance is calculated in the low-dimensional manifold space. The dominant direction of trajectory deviation is extracted based on the geodesic distance. The target tracking deviation is determined based on the geodesic distance and the dominant direction. The membership value corresponding to the target tracking deviation is calculated based on the preset Gaussian membership function. The target tracking deviation is mapped to a fuzzy linguistic variable based on the membership value. Fuzzy inference is performed on the fuzzy linguistic variable based on the preset fuzzy rule base to obtain the activation intensity of each fuzzy rule. Fuzzy control output is calculated based on the activation intensity and the fuzzy rule. The fuzzy control output is defuzzified using the centroid method to obtain the channel control component. The comprehensive driving command is obtained by superimposing the amplitude normalized components of the channel control components.

[0064] Feature analysis identifies and corrects the predicted trajectory to match the current task. Feature analysis includes extracting trajectory feature vectors and classifying them for the task. Each trajectory feature vector consists of ten features: the maximum, minimum, average, standard deviation, and number of peaks for the position trajectory, and five corresponding features for the force trajectory. For example, for a precision contact task, typical values ​​for the position trajectory feature vector are [50.2, 49.8, 50.0, 0.1, 1] mm, and for the force trajectory feature vector, they are [2.5, 0.0, 1.2, 0.8, 1] N. A support vector machine (SVM) classifier is used to classify the feature vectors, categorizing the current task into one of the preset task types. Common task types for the Z-axis module include precision contact, pressure detection, elasticity measurement, constant force control, and contour tracking. The classifier uses a radial basis function kernel with a kernel parameter of 0.1, achieving a classification accuracy of over 95%.

[0065] The target trajectory corresponding to the current task is extracted from a pre-defined trajectory library. The trajectory library stores standard trajectories for different task types, with multiple trajectory samples for each task type. The library is organized using a tree index, with the root node representing the task type and the leaf nodes representing the specific trajectory data. For example, a precision contact task stores 10 standard trajectories, each containing 50 position-force pairs at 50 time points. The selection of the target trajectory is based on similarity matching; the Euclidean distance between the current trajectory and similar trajectories in the library is calculated, and the trajectory with the smallest distance is selected as the target trajectory. If the Euclidean distance between the corrected predicted trajectory for a precision contact task and a sample in the trajectory library is 0.15, which is less than the pre-defined threshold of 0.2, then that sample is selected as the target trajectory.

[0066] The corrected predicted trajectory pairs and the target trajectory are embedded into a low-dimensional manifold space using a manifold learning algorithm. Manifold learning employs an isometric mapping algorithm to map high-dimensional trajectory data to a three-dimensional manifold space, preserving the geodesic distance relationships between data points. The algorithm parameters are set as follows: 20 nearest neighbors, 500 maximum iterations, and a convergence threshold of 0.001. For a position-force trajectory pair of length 50 (essentially a 100-dimensional vector), a three-dimensional representation is obtained after isometric mapping. For example, the corrected predicted trajectory for a precision contact mission has coordinates [0.3, 0.5, 0.2] in the manifold space, and the corresponding target trajectory has coordinates [0.4, 0.6, 0.3].

[0067] Geodesic distance is calculated in a low-dimensional manifold space. The geodesic distance calculation uses Dijkstra's algorithm, constructing a nearest neighbor graph on the manifold with edge weights equal to Euclidean distance, and solving for the shortest path length between two points. During the nearest neighbor graph construction, each point is connected to its 10 nearest neighbors. For the predicted trajectory and the target trajectory in the aforementioned example, the geodesic distance is calculated to be 0.18, representing the actual degree of deviation between the trajectories.

[0068] The dominant direction of trajectory deviation is extracted based on geodesic distance. The dominant direction is determined by calculating the gradient direction from the predicted trajectory point to the target trajectory. In the manifold space, the vector from the predicted point to the nearest point on the target trajectory is calculated, and this vector, after normalization, becomes the dominant direction. In the three-dimensional manifold space, the dominant direction is represented as a unit vector, for example, [0.707, 0.707, 0], indicating the dominance of deviation in the first two dimensions.

[0069] Target tracking deviation is determined based on geodesic distance and dominant direction. Target tracking deviation is the projection of geodesic distance onto the dominant direction and is calculated as the product of the geodesic distance and the components of the dominant direction. For the Z-axis module, the tracking deviation includes two components: position deviation and force deviation. Position deviation represents the difference between the current predicted position and the target position, for example, 0.12 mm; force deviation represents the difference between the current predicted force and the target force, for example, 0.35 Newtons.

[0070] The membership values ​​corresponding to the target tracking deviation are calculated based on a preset Gaussian membership function. The Gaussian membership function is determined by the center point and width parameters. For position deviation, five fuzzy sets are set: negative large (center -0.5 mm, width 0.2 mm), negative small (center -0.1 mm, width 0.1 mm), zero (center 0 mm, width 0.05 mm), positive small (center 0.1 mm, width 0.1 mm), and positive large (center 0.5 mm, width 0.2 mm). For force deviation, five fuzzy sets are similarly set, with parameters adjusted accordingly. The membership values ​​of position deviation and force deviation with respect to each fuzzy set are calculated. For example, a position deviation of 0.12 mm has a membership value of 0.85 with "positive small" and a membership value of 0.12 with "zero".

[0071] The target tracking deviation is mapped to a fuzzy linguistic variable based on the membership degree value. The fuzzy linguistic variable is the linguistic description corresponding to the fuzzy set with the largest membership degree. For example, the fuzzy linguistic variable for a position deviation of 0.12 mm is "positive small", and the fuzzy linguistic variable for a force deviation of 0.35 N is also "positive small".

[0072] Fuzzy inference calculations are performed on fuzzy linguistic variables based on a pre-defined fuzzy rule base. Fuzzy rules adopt the IF-THEN form, for example, "IF position deviation is positive and small AND force deviation is positive and small THEN position control is negative and small AND force control is negative and small". The fuzzy rule base for the Z-axis module contains 25 rules, covering various combinations of position and force deviations. The activation strength of a fuzzy rule is calculated as the minimum membership value of the antecedent part; for example, the activation strength of the aforementioned rule is min(0.85, 0.75) = 0.75.

[0073] The fuzzy control output is calculated based on activation intensity and fuzzy rules. Each activated fuzzy rule generates a control output suggestion, with the output intensity equal to the rule's activation intensity. For example, if a rule suggests position control as "negative small" and has an activation intensity of 0.75, then the rule's contribution to position control is the membership function "negative small" multiplied by 0.75. The contributions of all rules are summed to form the fuzzy control output, which is represented by two fuzzy sets: position control and force control.

[0074] The fuzzy control output is defuzzified using the centroid method. The centroid method calculates the centroid of the fuzzy set as the precise control value. Specifically, it is calculated by dividing the integral of the product of the fuzzy set and the membership value by the integral of the membership value. For example, the fuzzy output of position control becomes -0.08 mm after defuzzification, and force control becomes -0.25 N, representing the adjustment amount of position and force, respectively.

[0075] The integrated drive command is obtained by superimposing the channel control components after amplitude normalization. The normalization operation maps the control components to the range [-1, 1], dividing the position control component by the maximum permissible position adjustment (e.g., 0.5 mm) and the force control component by the maximum permissible force adjustment (e.g., 2 Newtons). The normalized position control is -0.16, and the force control is -0.125. The integrated drive command is obtained by weighted superposition of the two control components, with the weights set according to the task type. For precision contact tasks, the position control weight is 0.6, and the force control weight is 0.4, resulting in an integrated drive command of -0.16 × 0.6 + (-0.125) × 0.4 = -0.146.

[0076] In this embodiment, the current task corresponding to the corrected predicted trajectory is identified by feature analysis and matched with the target trajectory in the trajectory library. The reference trajectory can be automatically selected according to the task characteristics, which improves the targeting and task adaptability of trajectory tracking. The corrected predicted trajectory and the target trajectory are embedded into a low-dimensional manifold space by manifold learning algorithm and the deviation direction is calculated based on geodesic distance. This can more accurately reflect the real geometric difference of the trajectory in nonlinear space, improve the accuracy and robustness of trajectory deviation identification. Based on fuzzy rule base reasoning and combined with activation intensity weighting, flexible control output under multi-rule parallel response is realized, which effectively avoids the response rigidity problem caused by fixed parameters.

[0077] In one alternative implementation, The Z-axis module is driven according to the integrated drive command, and the stability index of the drive response is monitored. If the stability index exceeds a preset stability threshold, the amplitude and rate of change of the integrated drive command are limited to obtain a safe drive command and executed, including: The integrated drive command is converted into a pulse width modulation signal and applied to the drive motor of the Z-axis module. Real-time force response signal and real-time position response signal are collected. The real-time force response signal and real-time position response signal are decomposed into multiple layers by wavelet packet decomposition algorithm to extract energy value and construct energy feature vector. The center of the minimum enclosing hypersphere of the energy feature vector is calculated by support vector data description algorithm. The Mahalanobis distance from the energy feature vector to the center of the sphere is calculated to obtain the stability index. When the stability index is greater than the preset stability threshold, the amplitude of the integrated drive command is extracted and compared with the preset amplitude upper limit. The amplitude exceeding the preset amplitude upper limit is clipped to the preset amplitude upper limit to obtain the amplitude clipped command. The difference between the amplitude clipped command at adjacent time points is calculated to obtain the amplitude change rate. The amplitude change rate is compared with the preset change rate upper limit and the part exceeding the preset change rate upper limit is smoothed to obtain the safety drive command. The safety drive command is converted into a pulse width modulation signal and applied to the drive motor.

[0078] The integrated drive command is converted into a pulse width modulation (PWM) signal and applied to the drive motor of the Z-axis module. The integrated drive command value ranges from [-1, 1], and is linearly mapped to the duty cycle of the PWM signal. The mapping relationship is: duty cycle = (integrated drive command + 1) × 50%. For example, when the drive command is 0.5, the corresponding duty cycle is 75%; when the drive command is -0.8, the corresponding duty cycle is 10%. The frequency of the PWM signal is set to 20 kHz to ensure smooth motor response and moderate heat generation. The PWM signal is generated by the timer module of the microcontroller and output to the power drive circuit to drive the Z-axis module motor. The motor drive circuit adopts an H-bridge structure with an input voltage of 24 volts and a maximum current of 5 amps, which can meet the power requirements of the Z-axis module under various operating conditions.

[0079] Real-time force and position response signals are acquired. The force response signal is acquired by a force sensor installed at the end of the Z-axis module, with a sampling frequency of 1000 Hz, a measurement range of 0-100 Newtons, and a resolution of 0.01 Newtons. The position response signal is acquired by an optical encoder integrated on the Z-axis module, with an encoder accuracy of 0.001 mm and a sampling frequency of 1000 Hz. The acquired signals are pre-amplified and anti-aliasing filtered before being converted into digital signals by an analog-to-digital converter. The latency of the data acquisition system is controlled to within 1 millisecond to ensure timely response to changes in status.

[0080] The real-time force response signal and real-time position response signal are decomposed into energy values ​​through multi-level decomposition using a wavelet packet decomposition algorithm to construct an energy feature vector. The wavelet packet decomposition employs the Daubechies4 wavelet basis function, performing a four-level decomposition to generate 16 frequency band sub-signals. Each frequency band sub-signal contains the system's dynamic characteristics across different frequency ranges. For example, the first level of decomposition divides the 0-500 Hz signal into two parts: 0-250 Hz and 250-500 Hz; the second level further decomposes it into four parts: 0-125 Hz, 125-250 Hz, 250-375 Hz, and 375-500 Hz, and so on. The energy value, i.e., the sum of squares of the signal, is calculated for each frequency band sub-signal. The 16 frequency band energy values ​​of the force response signal and the 16 frequency band energy values ​​of the position response signal are combined to form a 32-dimensional energy feature vector. For example, the energy values ​​of the force response signal collected at a certain moment after decomposition are [10.5, 8.2, 5.1, 4.8, 3.2, 2.5, 1.8, 1.5, 1.2, 1.0, 0.8, 0.6, 0.4, 0.3, 0.2, 0.1], and the energy values ​​of the position response signal are [12.3, 9.5, 6.2, 4.9, 3.5, 2.8, 2.0, 1.7, 1.4, 1.1, 0.9, 0.7, 0.5, 0.4, 0.3, 0.2]. These values ​​are combined to form an energy feature vector.

[0081] The center of the minimum enclosing hypersphere of energy feature vectors is calculated using the Support Vector Data Description (SDR) algorithm. SDR is an algorithm used for anomaly detection by finding the minimum hypersphere that contains most of the normal data. The algorithm parameters are set as follows: a Gaussian radial basis function kernel with a kernel parameter of 0.5 and a relaxation parameter of 0.05. During the training phase, 1000 energy feature vector samples from the Z-axis module under normal operating conditions are collected to train the SDR model, obtaining the center and radius of the enclosing hypersphere. The center is a 32-dimensional vector representing the central location of the energy distribution under normal operating conditions. For example, the first 8 dimensions of the trained center are [11.2, 8.8, 5.5, 4.8, 3.4, 2.6, 1.9, 1.6]. The radius is typically set so that 95% of the normal samples are located inside the hypersphere, for example, a radius of 3.2.

[0082] The stability index is obtained by calculating the Mahalanobis distance from the energy eigenvector to the center of the sphere. Mahalanobis distance considers the correlation between features and reflects the data distribution in multidimensional space better than Euclidean distance. The calculation method involves multiplying the difference vector between the energy eigenvector and the center of the sphere by the inverse of the covariance matrix, then multiplying by the transpose of the difference vector, and finally taking the square root. The covariance matrix is ​​calculated using sample data from normal operating conditions. For example, a Mahalanobis distance of 2.8 from the energy eigenvector to the center of the sphere at a certain moment indicates the degree of deviation between the current state and the normal state. The stability index is defined as this Mahalanobis distance; the smaller the value, the more stable the system.

[0083] When the stability index exceeds a preset stability threshold, a security control mechanism is triggered. The stability threshold is set based on the application scenario, with a typical value of 3.5. This value is determined through statistical analysis of normal and abnormal states, ensuring that 99% of normal states will not trigger security controls and that potential unstable states can be detected promptly. For example, if the current stability index is 3.8, which is greater than the threshold of 3.5, it is determined that the current state may be unstable and security control measures are required.

[0084] Extract the amplitude of the integrated drive command and compare it with the preset amplitude upper limit. The amplitude equals the absolute value of the integrated drive command. The preset amplitude upper limit is determined based on the safe operating range of the Z-axis module and is usually set to 0.8. If the current integrated drive command is -0.95, its amplitude of 0.95 is greater than the preset upper limit of 0.8, amplitude clipping is required.

[0085] The amplitude exceeding the preset amplitude limit is clipped to the preset amplitude limit to obtain the clipped amplitude command. The clipping operation keeps the sign of the command unchanged, only adjusting the amplitude. For example, the original integrated drive command is -0.95, and the clipped value is -0.8; the original command is 0.9, and the clipped value is 0.8. Amplitude clipping ensures that the drive command does not exceed the safe operating range of the Z-axis module, preventing motor overload or mechanical damage.

[0086] The amplitude change rate is obtained by calculating the difference between the commands after amplitude clipping at adjacent time points. The rate of change is the difference between the commands of two consecutive control cycles divided by the control cycle time. For example, if the clipped commands for two adjacent control cycles are 0.6 and 0.75 respectively, and the control cycle is 1 millisecond, then the rate of change is 150 units per second. An excessively high rate of change can cause violent movement of the Z-axis module, leading to mechanical vibration or even damage.

[0087] The amplitude change rate is compared with a preset upper limit, and the portion exceeding the upper limit is smoothed. The preset upper limit is typically set to 100 units per second. The smoothing process uses a sliding window averaging method with a window size of 5 sampling points. For example, if the current change rate is 150 units per second, exceeding the upper limit of 100 units per second, the five most recent control commands are weighted and averaged to obtain a smoothed command, reducing the change rate to a safe range. The smoothed command provides a safe drive command, ensuring smooth and safe movement of the Z-axis module.

[0088] The safety drive command is converted into a pulse width modulation signal and applied to the drive motor, and the safety drive command is linearly mapped to a duty cycle. For example, a safety drive command of 0.65 corresponds to a duty cycle of 82.5%.

[0089] In this embodiment, the response signal is decomposed into multiple energy layers using a wavelet packet decomposition algorithm, which can simultaneously capture subtle fluctuations and low-frequency trends across different frequency bands, constructing a more comprehensive energy feature representation. The minimum enclosing hypersphere of the energy feature vector is calculated using a support vector data description algorithm, and a stability index is defined based on Mahalanobis distance. This allows for a quantitative assessment of the deviation from the current operating state, providing early warning of abnormal vibrations or dynamic instability. The drive command is adaptively corrected through a dual constraint mechanism of amplitude clipping and rate of change smoothing, effectively suppressing abrupt changes and overshoot in the drive signal, preventing structural impacts and control oscillations caused by motor overload or severe acceleration / deceleration. By converting the smoothed and corrected safety drive command back into a pulse width modulation signal and applying it to the drive motor, the entire control closed loop possesses dynamic stability protection, significantly improving the operational stability and safety of the Z-axis module under complex dynamic loads, reducing mechanical fatigue of the actuator and energy consumption fluctuations in the control system, and achieving a synergistic unity of high-precision drive and safety constraints.

[0090] A second aspect of the present invention provides a high-precision Z-axis module integrated drive control system based on force feedback, comprising: The first unit is used to acquire the position signal and force sensing signal of the Z-axis module during the motion process and synchronize them according to the time sequence, and calculate a multi-dimensional state sequence based on the position signal and the force sensing signal. The second unit is used to calculate the degree of correlation between the components of each dimension in the multidimensional state sequence to obtain the coupling relationship and determine the coupling influence coefficient. Based on the coupling influence coefficient, the multidimensional state sequence is linearly transformed to obtain the channel decoupled state sequence. The channel decoupled state sequence is segmented by a sliding window and the local change trend is calculated. Based on the local change trend, the position evolution trajectory and the force evolution trajectory are predicted and time-series aligned to obtain the predicted state trajectory pair. The third unit is used to perform cross-correlation analysis on the predicted state trajectory pair, identify the bidirectional coupling influence relationship between position and force and construct a bidirectional causal chain, dynamically adjust the predicted state trajectory pair based on the bidirectional causal chain to obtain a corrected predicted trajectory pair, calculate the target tracking deviation of the corrected predicted trajectory pair and solve for the channel control component, and superimpose the channel control component to obtain a comprehensive driving command. The fourth unit is used to apply drive to the Z-axis module according to the comprehensive drive command and monitor the stability index of the drive response. If the stability index exceeds the preset stability threshold, the comprehensive drive command is subjected to amplitude limitation and rate of change constraint to obtain a safe drive command and execute it.

[0091] A third aspect of the present invention provides an electronic device, comprising: A processor and a memory for storing processor-executable instructions, wherein the processor is configured to invoke instructions stored in the memory to perform the aforementioned method.

[0092] A fourth aspect of the present invention provides a computer-readable storage medium having stored thereon computer program instructions that, when executed by a processor, implement the aforementioned method.

[0093] This invention can be a method, apparatus, system, and / or computer program product. The computer program product may include a computer-readable storage medium having computer-readable program instructions loaded thereon for performing various aspects of the invention.

[0094] Finally, it should be noted that the above embodiments are only used to illustrate the technical solutions of the present invention, and not to limit them; although the present invention has been described in detail with reference to the foregoing embodiments, those skilled in the art should understand that modifications can still be made to the technical solutions described in the foregoing embodiments, or equivalent substitutions can be made to some or all of the technical features; and these modifications or substitutions do not cause the essence of the corresponding technical solutions to deviate from the scope of the technical solutions of the embodiments of the present invention.

Claims

1. A high-precision Z-axis module integrated drive control method based on force feedback, characterized in that, The method comprises the following steps: acquiring position signals and force sensing signals of a Z-axis module during movement and synchronously aligning them according to time sequences, calculating a multi-dimensional state sequence based on the position signals and the force sensing signals; calculating the coupling relationship between the dimension components in the multi-dimensional state sequence to obtain a coupling influence coefficient, performing linear transformation on the multi-dimensional state sequence based on the coupling influence coefficient to obtain a channel decoupling state sequence, performing sliding window segmentation on the channel decoupling state sequence and calculating a local change trend, predicting a position evolution trajectory and a force evolution trajectory based on the local change trend and performing time sequence alignment to obtain a predicted state trajectory pair; performing cross-correlation analysis on the predicted state trajectory pair, identifying the bidirectional coupling influence relationship between the position and the force and constructing a bidirectional causal chain, performing dynamic constraint adjustment on the predicted state trajectory pair based on the bidirectional causal chain to obtain a modified predicted trajectory pair, calculating a target tracking deviation of the modified predicted trajectory pair and solving to obtain a channel control component, and superimposing the channel control component to obtain a comprehensive driving instruction; applying driving to the Z-axis module according to the comprehensive driving instruction and monitoring a stability index of the driving response, if the stability index exceeds a preset stability threshold, performing amplitude limitation and change rate constraint on the comprehensive driving instruction to obtain a safe driving instruction and executing the safe driving instruction.

2. The method of claim 1, wherein, The method comprises the following steps: acquiring position signals and force sensing signals of a Z-axis module during movement and synchronously aligning them according to time sequences, calculating a multi-dimensional state sequence based on the position signals and the force sensing signals; collecting the position signals of the Z-axis module through a position sensor and collecting the force sensing signals of the Z-axis module through a force sensor, extracting the timestamp information of the position signals and the force sensing signals respectively and identifying the time offset between the position signals and the force sensing signals, performing time sequence interpolation alignment based on the time offset to obtain aligned position signals and aligned force signals; performing first-order differential operation on the aligned position signals to obtain instantaneous velocity signals, performing second-order differential operation on the aligned position signals to obtain instantaneous acceleration signals, combining the instantaneous velocity signals and the instantaneous acceleration signals into motion derivative components, and performing high-frequency noise suppression processing on the aligned force signals to extract effective frequency band components as net force components; 3. The method of claim 1, wherein, vector splicing the aligned position signals, the motion derivative components and the net force components according to the time dimension to obtain the multi-dimensional state sequence. The method comprises the following steps: performing cross-correlation operation on any two dimension components in the multi-dimensional state sequence to obtain a correlation coefficient matrix, performing eigenvalue decomposition on the correlation coefficient matrix, extracting a feature vector with an eigenvalue greater than a preset coupling judgment threshold as a target feature vector, and identifying a dimension combination with a coupling relationship based on the target feature vector to obtain a coupling dimension combination. construct a joint state space based on the position dimension component and the force dimension component in the coupling dimension combination, perform phase space reconstruction on the joint state space to obtain a coupling phase space trajectory, and perform dynamic time warping analysis on the coupling phase space trajectory to extract a time alignment mapping relationship between an evolution path of the position dimension component and an evolution path of the force dimension component; perform bidirectional regression fitting on the position dimension component and the force dimension component based on the time alignment mapping relationship to obtain a pair of bidirectional transfer functions, wherein the pair of bidirectional transfer functions constitute a coupling relationship, and perform frequency domain analysis on the pair of bidirectional transfer functions to obtain a frequency response characteristic and extract a gain coefficient of a dominant frequency component as a coupling influence coefficient; construct a decoupling transformation matrix based on the coupling influence coefficient and a target feature vector, and perform matrix left multiplication transformation on a multi-dimensional state sequence to obtain a position channel state sequence and a force channel state sequence and combine the position channel state sequence and the force channel state sequence to obtain a channel decoupled state sequence.

4. The method of claim 1, wherein, perform sliding window segmentation on the channel decoupled state sequence and calculate a local change trend, predict a position evolution trajectory and a force evolution trajectory based on the local change trend, and perform time sequence alignment to obtain a predicted state trajectory pair including: set a sliding window for the channel decoupled state sequence, perform sliding segmentation on the channel decoupled state sequence according to a preset step size to obtain a plurality of window segments, perform multi-dimensional linear fitting on data in each window segment, calculate a fitting slope of a position channel dimension to obtain a position change rate, calculate a fitting slope of a force channel dimension to obtain a force change rate, and combine the position change rate and the force change rate to obtain a local change trend; based on the local change trend, extract an end value of the position channel in the channel decoupled state sequence as a position starting point, linearly extrapolate from the position starting point at an extension rate of the position change rate to obtain a position evolution trajectory, extract an end value of the force channel in the channel decoupled state sequence as a force starting point, and linearly extrapolate from the force starting point at an extension rate of the force change rate to obtain a force evolution trajectory; extract a first timestamp sequence corresponding to the position evolution trajectory and a second timestamp sequence corresponding to the force evolution trajectory, perform time sequence alignment on the first timestamp sequence and the second timestamp sequence to obtain a unified time reference, map the position evolution trajectory and the force evolution trajectory to the same time coordinate system based on the unified time reference, and obtain a predicted state trajectory pair.

5. The method of claim 1, wherein, perform cross-correlation analysis on the predicted state trajectory pair, identify a bidirectional coupling influence relationship between the position and the force, and construct a bidirectional causal chain, perform dynamic constraint adjustment on the predicted state trajectory pair based on the bidirectional causal chain, and obtain a modified predicted trajectory pair including: extract a position value sequence corresponding to the position evolution trajectory pair and a force value sequence corresponding to the force evolution trajectory pair from the predicted state trajectory pair, perform a differential operation on the position value sequence to obtain a velocity sequence, perform a differential operation on the velocity sequence to obtain an acceleration sequence, calculate a cross-correlation function of the acceleration sequence and the force value sequence, extract an inverse of a dynamic response delay of force to position at a peak time of the cross-correlation function to obtain a force causal strength, perform an integral operation on the force value sequence to obtain an impulse sequence, calculate a partial mutual information of the impulse sequence and the position value sequence and extract an energy transfer efficiency of position to force to obtain a position causal strength, and identify a bidirectional coupling influence relationship based on the force causal strength and the position causal strength; set a position state and a force state as nodes of a causal graph, establish a directed connection from a position node to a force node and a directed connection from the force node to the position node based on the bidirectional coupling influence relationship, and combine the directed connections to obtain a bidirectional causal chain; perform path analysis on the bidirectional causal chain to identify a forward propagation path from the position node to the force node and a reverse propagation path from the force node to the position node and establish a bidirectional coupling constraint relationship, and correct the predicted state trajectory pair based on the bidirectional coupling constraint relationship to obtain a corrected predicted trajectory pair.

6. The method of claim 1, wherein, calculate a target tracking deviation of the corrected predicted trajectory pair and solve to obtain a channel control component, and superimpose the channel control component to obtain a comprehensive driving instruction including: identify a current task corresponding to the corrected predicted trajectory pair through feature analysis, extract a target trajectory corresponding to the current task from a preset trajectory library, embed the corrected predicted trajectory pair and the target trajectory into a low-dimensional manifold space through a manifold learning algorithm and calculate a geodesic distance in the low-dimensional manifold space, extract a dominant direction of trajectory deviation based on the geodesic distance, and determine a target tracking deviation based on the geodesic distance and the dominant direction; calculate a membership value corresponding to the target tracking deviation based on a preset Gaussian membership function, map the target tracking deviation to a fuzzy language variable based on the membership value, calculate an activation strength of each fuzzy rule based on a preset fuzzy rule base through fuzzy reasoning of the fuzzy language variable, calculate a fuzzy control output based on the activation strength and the fuzzy rule, and de-fuzzify the fuzzy control output through a barycenter method to obtain a channel control component; perform amplitude normalization on the channel control component and superimpose to obtain the comprehensive driving instruction.

7. The method of claim 1, wherein, apply driving to the Z-axis module according to the comprehensive driving instruction and monitor a stability index of a driving response, if the stability index exceeds a preset stability threshold, perform amplitude limitation and a change rate constraint on the comprehensive driving instruction to obtain a safe driving instruction and execute including: The integrated driving instruction is converted into a pulse width modulation signal and applied to a driving motor of a Z-axis module, real-time force response signals and real-time position response signals are collected, energy values are extracted by multi-layer decomposition of the real-time force response signals and the real-time position response signals through a wavelet packet decomposition algorithm to construct an energy feature vector, a minimum enclosing hypersphere of the energy feature vector is calculated through a support vector data description algorithm, and a Mahalanobis distance of the energy feature vector to the center of the hypersphere is calculated to obtain a stability index; When the stability index is greater than a preset stability threshold, the amplitude of the integrated driving instruction is extracted and compared with a preset amplitude upper limit, the amplitude exceeding the preset amplitude upper limit is clipped to the preset amplitude upper limit to obtain an amplitude-clipped instruction, a differential value of the amplitude-clipped instruction at adjacent time instants is calculated to obtain an amplitude change rate, the amplitude change rate is compared with a preset change rate upper limit, and a portion exceeding the preset change rate upper limit is smoothed to obtain a safe driving instruction, and the safe driving instruction is converted into a pulse width modulation signal and applied to the driving motor.

8. A high-precision Z-axis module integrated drive control system based on force feedback, for implementing the method of any one of the preceding claims 1-7, characterized in that, Comprise: A first unit for acquiring position signals and force sensing signals of a Z-axis module during movement and synchronously aligning in time sequence, and calculating a multi-dimensional state sequence based on the position signals and the force sensing signals; A second unit for calculating a coupling relationship between the dimensions in the multi-dimensional state sequence to obtain a coupling influence coefficient, performing linear transformation on the multi-dimensional state sequence based on the coupling influence coefficient to obtain a channel decoupling state sequence, segmenting the channel decoupling state sequence by a sliding window and calculating a local change trend, predicting a position evolution trajectory and a force evolution trajectory based on the local change trend and performing time sequence alignment to obtain a predicted state trajectory pair; A third unit for cross-correlation analysis of the predicted state trajectory pair, identifying a bidirectional coupling influence relationship between the position and the force and constructing a bidirectional causal chain, performing dynamic constraint adjustment on the predicted state trajectory pair based on the bidirectional causal chain to obtain a modified predicted trajectory pair, calculating a target tracking deviation of the modified predicted trajectory pair and solving to obtain a channel control component, and superimposing the channel control component to obtain an integrated driving instruction; A fourth unit for applying driving to the Z-axis module according to the integrated driving instruction and monitoring a stability index of the driving response, and performing amplitude limitation and change rate constraint on the integrated driving instruction if the stability index exceeds a preset stability threshold to obtain a safe driving instruction and execute the safe driving instruction.

9. An electronic device, comprising: Comprise: A processor; A memory for storing processor-executable instructions; The processor is configured to call the instructions stored in the memory to execute the method of any one of claims 1 to 7.

10. A computer-readable storage medium having stored thereon computer program instructions, wherein, The computer program instructions are executed by the processor to implement the method of any one of claims 1 to 7.

Images