Industrial stacker intelligent control system based on machine vision

Abstract

The invention relates to the technical field of industrial control, in particular to an industrial stacker intelligent control system based on machine vision, and the system comprises a state resolving module which is used for carrying out the feature region decoupling of original observation data, and resolving a position deviation matrix and a deformation state variable for a rigid subspace and a flexible subspace respectively; the dynamic modeling module is used for identifying the first-order inherent frequency and the damping ratio of the mechanical structure on line and establishing a second-order impedance contact model; the impedance control module is used for generating an expected speed based on the position deviation and adaptively adjusting an impedance parameter according to a deformation state so as to carry out compliant modulation on a speed vector; and the estimation shaping module is used for correcting system time lag by using a Kalman observer and inhibiting residual vibration in combination with a zero-vibration input shaper. Through multi-dimensional physical field fusion and dynamic modeling, the problem that a rigid control model is mismatched with a flexible operation object is solved, and self-adaptive compliant control and vibration-free accurate operation of the stacking machine in the unstructured environment are achieved.

Description

Technical Field

[0001] This invention relates to the field of industrial control technology, specifically to an intelligent control system for industrial stacker cranes based on machine vision. Background Technology

[0002] In the fields of industrial process control and motion control, servo systems for long-stroke precise positioning of large inertia loads typically employ a semi-closed-loop position control architecture based on motor encoder feedback. However, such control systems suffer from the following technical drawbacks when dealing with complex controlled objects exhibiting time-varying parameters or unstructured disturbances: In the field of industrial motion control, existing servo control systems are mostly designed based on linear time-invariant theory. Their constant control parameters after initialization cannot adapt to the nonlinear drift of load inertia and damping ratio during operation. This model mismatch, due to the lack of online adaptive adjustment mechanism, leads to system overshoot or steady-state error, which seriously undermines the robustness of the closed-loop system. At the same time, traditional rigid body motion planning algorithms ignore the unmodeled dynamics introduced by the structural elasticity of large motion systems, and cannot effectively suppress low-frequency residual vibrations excited by inertial forces, forcing the system to maintain stability at the expense of bandwidth. In addition, when constructing full closed-loop control, the multi-rate sampling characteristics between the external observation unit and the underlying servo controller introduce nondeterministic pure time delay elements into the feedback loop, directly reducing the phase margin of the system. If time delay data is directly introduced for feedback correction, it is very easy to cause the closed-loop poles to shift out of the unit circle, thereby triggering continuous limit cycle oscillations.

[0003] To address the technical challenge of uncontrolled contact force caused by time-varying deviations between end effectors and flexible work objects in highly dynamic, unstructured warehousing environments, a machine vision-based intelligent control system for industrial stacker cranes is proposed. Summary of the Invention

[0004] The purpose of this invention is to provide an intelligent control system for industrial stacker cranes based on machine vision. It aims to achieve compliant adaptive control and high-dynamic vibration-free precision operation under nonlinear time-varying load conditions by integrating multi-dimensional state observation, adaptive impedance adjustment and dynamic shaping control technology.

[0005] An adaptive control system for an industrial stacker crane based on multi-dimensional state observation includes: State calculation module: Performs feature region decoupling processing on the raw observation data, performs 6D pose calculation for the rigid state subspace to obtain the position deviation matrix, and performs manifold geometry calculation for the flexible state subspace to obtain the deformation state variables; Dynamic modeling module: Monitors the displacement response sequence of the top of the stacker crane column, identifies the first-order natural frequency and damping ratio online; establishes a second-order impedance contact model including virtual mass, damping coefficient and stiffness coefficient; Impedance control module: maps the position deviation matrix to the desired velocity vector of the servo system; adaptively adjusts the stiffness coefficient and damping coefficient of the second-order impedance contact model according to the deformation state variables, and uses the adjusted impedance contact model to perform compliance modulation on the desired velocity vector to generate initial control commands. The estimation and shaping module constructs a Kalman state observer to perform state estimation on the initial control command and generate a corrected velocity control command. Based on the identified first-order natural frequency and damping ratio, a zero-vibration input shaper is constructed to perform time-domain convolution modulation on the corrected velocity control command and generate a shaping execution command.

[0006] Preferably, the acquisition of the raw observation data includes: A time-division trigger signal is sent to the dual-channel optical sensing unit. First, the first channel is activated to project near-infrared speckle structured light, and the depth field information is analyzed using the long-wavelength low scattering characteristics of the near-infrared spectrum. Then, the second channel is activated to project linearly polarized blue light, and in conjunction with the orthogonal analyzer at the imaging end, the surface optical property information is acquired using the principle of polarization extinction. Based on the pre-calibrated intrinsic parameter matrix and hand-eye extrinsic parameter matrix, the depth field information is converted into a three-dimensional spatial coordinate model, and the distortion of the surface optical property information is corrected. Spatial coordinate registration is performed to map the surface optical property values ​​to the corresponding three-dimensional spatial coordinate points, generating spatial state data as the original observation data.

[0007] Preferably, the state calculation module is specifically configured as follows: The three-dimensional spatial coordinates and surface optical attribute values ​​of the original observation data are extracted, a feature tensor is established, and input into the semantic segmentation neural network model to extract local topological features and surface attribute features. The state category probability of each data point is calculated. According to the maximum a posteriori probability principle, a state label is assigned to each data point to generate rigid and flexible classification masks. The flexible classification mask is applied to perform spatial filtering on the original observation data to separate the rigid state subspace and the flexible state subspace.

[0008] Preferably, the state calculation module is specifically configured to perform the calculation as follows: For the rigid state subspace, an iterative nearest-point registration algorithm is used to align the spatial data within the region with the pre-set rigid body standard model, and the rotation and translation transformation matrix of the current coordinate system relative to the ideal working coordinate system is calculated to obtain the position deviation matrix. For the flexible state subspace, an elliptical cylindrical manifold fitting algorithm is used to reconstruct the surface geometric topology of the region, calculate the ratio of major and minor axes and the principal axis tilt angle parameters of the flexible working object, and compare them with the pre-set object standard shape to obtain the deformation state variables.

[0009] Preferably, the dynamic modeling module is specifically configured as follows: A time-domain response sequence is constructed by capturing the horizontal vibration displacement at the top of the column. The first-order natural frequency is extracted by performing a fast Fourier transform, and the damping ratio is estimated using the logarithmic decay method. A second-order linear differential equation describing the dynamic mapping between the position correction and the contact force is established as the impedance contact model, in which the virtual mass, damping coefficient, and stiffness coefficient are set to characterize the inertial response, energy dissipation rate, and elastic restoring force of the system, respectively.

[0010] Preferably, the impedance control module is specifically configured as follows: A position-based visual servo control law is used to construct a Jacobian matrix. The position deviation matrix is ​​used as feedback input. The desired velocity vector of each servo axis is calculated through inverse kinematics calculation and proportional gain adjustment. A nonlinear attenuation mapping function between deformation state variables and stiffness coefficients is established. When the deformation state variables increase, the driving stiffness coefficient decreases along the mapping function curve. Based on the constant damping ratio constraint, the damping coefficient is synchronously adjusted according to the real-time changing stiffness coefficient. An admittance control structure is used. The desired velocity vector is used as a reference input. The adjusted stiffness coefficient and damping coefficient are substituted into the second-order impedance contact model to construct a dynamic filter. The desired velocity vector is filtered, and the position and velocity corrections are calculated to generate the initial control command.

[0011] Preferably, the estimation and shaping module is specifically configured as follows: A kinematic recursive model is constructed using feedback data from a high-frequency sampled servo encoder. A Kalman filter is then used to perform time update calculations to obtain prior state estimates. A low-frequency sampled position deviation matrix is ​​used as the observation vector, and a Kalman filter is used to perform measurement update calculations. The prior state estimates are corrected using visual observation residuals to generate a corrected speed control command. The pulse time interval and amplitude coefficient are calculated using the identified first-order natural frequency and damping ratio to construct a zero-vibration pulse sequence. This zero-vibration pulse sequence is then used as a kernel function to perform a time-domain convolution operation on the corrected speed control command, resulting in a shaped execution command containing a stepped pulse sequence.

[0012] Preferably, the impedance control module further includes virtual potential field guiding logic, specifically configured as follows: The surface geometry is received, and a virtual repulsive potential field is constructed. The maximum contour envelope of the flexible cargo is defined as the zero-distance boundary of the potential field. A positive correlation mapping between the potential field intensity gradient and the deformation state variable is established. When the deformation state variable increases, the repulsive intensity of the potential field is dynamically increased. During the generation of the initial control command, the Euclidean distance gradient between the end effector and the potential field boundary is calculated in real time. A virtual repulsive velocity vector is generated and superimposed on the desired velocity vector. When the end effector approaches the deformation boundary of the cargo, a virtual reverse damping force is generated.

[0013] Compared with the prior art, the beneficial effects of the present invention are as follows: 1. By decoupling the state of the original observation data, the system can simultaneously ensure the positional alignment accuracy for rigid pallets and the safety of contact force control for flexible goods. By adaptively adjusting the stiffness and damping coefficient of the second-order impedance model using deformation state variables, the control law can smoothly switch from rigid position control to compliant force-position hybrid control. This effectively solves the problem of physical damage caused by overshoot or compression when gripping easily deformable goods such as tires using traditional rigid control, and significantly improves the system's process adaptability to non-standard goods.

[0014] 2. Breaking through the limitations of traditional kinematic control, an online dynamic identification mechanism is introduced, enabling real-time capture of the frequency drift of the column under different load conditions. A zero-vibration input shaper, built based on the identification results, decomposes the original step command into a self-canceling pulse sequence using time-domain convolution modulation technology, actively canceling the elastic sway of the mechanical structure. This allows the stacker crane to quickly stabilize after stopping at high speed, significantly reducing waiting time and improving warehouse throughput efficiency while ensuring positioning accuracy.

[0015] 3. To address the time delay issue caused by multi-rate sampling in industrial settings, the system employs a Kalman state estimation strategy. By fusing servo high-frequency feedback with visual low-frequency observation, the Kalman observer effectively compensates for the phase lag caused by system time delay, providing optimal state estimation. Furthermore, by combining virtual potential field guidance logic, a non-contact safety boundary based on deformation state is constructed at the algorithm level, generating a virtual reverse damping force. This ensures the convergence of the control system under strong interference and communication delay environments, preventing the risk of physical collisions. Attached Figure Description

[0016] Figure 1 This is a system architecture diagram of an intelligent control system for industrial stacker cranes based on machine vision, according to the present invention. Figure 2 This is a flowchart illustrating the execution logic of the state calculation module in this invention. Figure 3 This is a flowchart illustrating the intelligent control logic of the industrial stacker crane of the present invention. Detailed Implementation

[0017] 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.

[0018] Please see Figures 1 to 3This invention provides an intelligent control system for industrial stacker cranes based on machine vision. The technical solution is as follows, please refer to... Figure 1 For the system architecture diagram, refer to Figure 3 This is a flowchart, specifically including: State calculation module: Performs feature region decoupling processing on the raw observation data, performs 6D pose calculation for the rigid state subspace to obtain the position deviation matrix, and performs manifold geometry calculation for the flexible state subspace to obtain the deformation state variables; Dynamic modeling module: Monitors the displacement response sequence of the top of the stacker crane column, identifies the first-order natural frequency and damping ratio online; establishes a second-order impedance contact model including virtual mass, damping coefficient and stiffness coefficient; Impedance control module: maps the position deviation matrix to the desired velocity vector of the servo system; adaptively adjusts the stiffness coefficient and damping coefficient of the second-order impedance contact model according to the deformation state variables, and uses the adjusted impedance contact model to perform compliance modulation on the desired velocity vector to generate initial control commands. The estimation and shaping module constructs a Kalman state observer to perform state estimation on the initial control command and generate a corrected velocity control command. Based on the identified first-order natural frequency and damping ratio, a zero-vibration input shaper is constructed to perform time-domain convolution modulation on the corrected velocity control command and generate a shaping execution command. Example

[0019] In this embodiment, the system described in this application is applied to a gripping scenario involving metal connecting rods or precision gears with high specular reflection characteristics in an industrial automation implementation environment. The specific steps are as follows: The acquisition and processing of raw observation data is performed by an intelligent sensing controller deployed on the robot's end effector. This controller interacts with the dual-channel optical sensing unit and the robotic arm control system via gigabit Ethernet. Specifically, the step of sending a time-division trigger signal to the dual-channel optical sensing unit is strictly executed according to a preset timing logic. The timing generator inside the controller generates a synchronous acquisition cycle with a period of 50ms. At the beginning of each acquisition cycle (T=0), the controller outputs a first high-level trigger pulse with a pulse width set to 100μs. This signal is transmitted to the drive circuit of the first optical channel through a shielded twisted pair cable.

[0020] Furthermore, after the first optical channel is activated, the projector emits near-infrared speckle structured light with a center wavelength of 850nm. This band is chosen to utilize the long-wavelength, low-scattering characteristics of near-infrared light propagating in air, thereby enhancing the penetration stability of the beam in the oil mist environment of the workshop. At this time, the near-infrared photosensitive chip at the imaging end simultaneously activates the exposure, with the exposure time set to 5ms. The imaging unit acquires the reflected image with the speckle pattern, and the image resolution is set to 1280×960 pixels. In order to analyze the depth field information, the processor adopts a stereo matching algorithm based on local regions, setting the matching window size to 15×15 pixels. Within the preset parallax search range, the processor determines the corresponding point by calculating the normalized cross-correlation coefficient (NCC) of the pixel grayscale values ​​of the left and right image blocks. Then, based on the principle of triangulation, the product of the baseline distance (e.g., set to 120mm) and the lens focal length is divided by the calculated parallax value to obtain the vertical depth value of each pixel, generating a depth field matrix.

[0021] Subsequently, after a 10ms interval following the completion of the first channel acquisition, the controller outputs a second trigger signal to activate the second channel. The second channel projects linearly polarized blue light with a center wavelength of 450nm. Blue light is chosen because its short wavelength allows for more detailed micro-texture feedback on metal surfaces. At this time, the polarization direction of the projected light is set to the vertical direction, while the transmission axis of the analyzer configured in front of the imaging lens is set to the horizontal direction, forming an orthogonal state (90-degree angle). Utilizing the principle of polarization extinction, the directly reflected strong specular light is physically blocked, and only the diffusely reflected light that undergoes depolarization can pass through the analyzer and be captured by the blue light sensor chip. This step filters out high-light noise through physical means, acquiring surface optical property information that reflects the true texture and material characteristics of the object's surface, i.e., a grayscale intensity image, with a data bit depth set to 8 bits.

[0022] Furthermore, after acquiring the depth field information and surface optical property information, the system calls the calibration data pre-stored in the non-volatile memory. Among them, the intrinsic parameter matrix includes the focal length, principal point coordinates, and three radial distortion coefficients and two tangential distortion coefficients; the hand-eye extrinsic parameter matrix describes the rotation and translation relationship of the camera coordinate system relative to the robot end flange coordinate system; for the surface optical property information (i.e., polarized blue light image), the processor uses bilinear interpolation to perform distortion correction; specifically, for each integer coordinate pixel in the corrected image, the processor uses the distortion model to inversely calculate its floating-point coordinates in the original distorted image, and selects the gray values ​​of the four neighboring pixels around the floating-point coordinates, and performs a weighted sum according to the inverse proportional relationship of the distance weight from the pixel to the floating-point coordinate to obtain the corrected gray value.

[0023] Specifically, the process of performing spatial coordinate registration is to establish a unified data dimension. The processor back-projects each depth value in the depth field information, combined with the inverse matrix of the intrinsic parameter matrix, onto three-dimensional space to generate a three-dimensional point cloud in the camera coordinate system. Then, using the rotation matrix (3 rows and 3 columns) and translation vector (3 rows and 1 column) in the hand-eye extrinsic parameter matrix, the point cloud data is transformed into the robot base coordinate system through matrix multiplication and vector addition. At this point, since the depth sensor and polarization image sensor have undergone strict physical alignment and calibration, the system directly maps the corrected surface optical attribute information (grayscale values) onto the corresponding three-dimensional spatial coordinate points. If there is a physical positional deviation between the two sensors, the three-dimensional points are projected onto the polarization image plane according to the rigid body transformation matrix between the two sensors, and the grayscale value of the nearest neighbor pixel is selected as the attribute value. The final generated spatial state data is a four-dimensional tensor, containing the X-axis coordinate, Y-axis coordinate, Z-axis coordinate, and surface texture grayscale value of each point. Furthermore, after acquiring the aforementioned spatial state data, in order to achieve the identification and localization of the target object, this embodiment employs a feature parsing network based on deep learning to process the spatial state data; Specifically, the spatial state data is first structured into an input tensor of size 256×256×4; the first three channels correspond to the normalized values ​​of the X, Y, and Z spatial coordinates, and the fourth channel corresponds to the grayscale value of the polarization texture; this input tensor is fed into a custom convolutional neural network model, the specific configuration of which is as follows: The first layer is the feature encoding layer, which contains 32 convolutional kernels. The size of each convolutional kernel is set to 3×3 and the stride is set to 1. This layer uses the Modified Linear Unit (ReLU) as the activation function. In this step of the calculation, the convolutional kernel slides on the input tensor by the stride, performs a dot product operation on the spatial coordinates and texture grayscale within each 3×3 local window, and sums them up to extract the primary geometric edge features.

[0024] Subsequently, the data flows through a max pooling layer with a pooling window size of 2×2 and a stride of 2, which downsamples the spatial dimension of the feature map to 128×128 to reduce computation and extract translation-invariant features. Furthermore, the network includes a core feature extraction module consisting of three consecutive convolutional layers with the number of convolutional kernels increasing sequentially to 64, 128, and 256, while maintaining a kernel size of 3×3. After each convolutional operation, batch normalization is performed, adjusting the distribution of feature values ​​to a standard normal distribution with a mean of 0 and a variance of 1 to accelerate model convergence and prevent gradient vanishing. The model's output is connected to a fully connected layer containing 512 neurons, and overfitting is prevented using Dropout technology (with a dropout rate set to 0.5). Finally, a regression layer with an output dimension of 6 outputs the target's pose parameters, including three position coordinate values ​​and three Euler angle values.

[0025] In terms of processing logic, the system sets a confidence threshold of 0.85. When the confidence score of the pose parameters output by the network is higher than 0.85, the pose data is deemed valid and sent to the robot controller via the industrial bus to perform the grasping action. If it is lower than the threshold, the system will trigger an alarm and instruct the robot to readjust the observation angle, repeating the above time-division triggering and data acquisition process. Through this closed-loop data flow, the feasibility and reliability of robot operation in complex optical environments are ensured.

[0026] This application effectively solves the problem of 3D reconstruction on highly reflective metal surfaces or in oil mist environments in industrial settings by fusing near-infrared speckle structured light with linearly polarized blue light. It utilizes the principle of polarization extinction to physically filter out strong light noise caused by specular reflection, and combines the long-wavelength low-scattering penetration characteristics of near-infrared light to significantly improve the signal-to-noise ratio, integrity, and texture clarity of the original observation data under complex optical conditions.

[0027] Furthermore, the construction and preprocessing of the feature tensor are as follows: First, the processor performs structured sampling on the raw observation data generated in the previous steps to construct a feature tensor that meets the input requirements of the neural network. Since the number of raw observation data points is huge and unevenly distributed, the processor executes the farthest point sampling algorithm to select 4096 data points that are representative and spatially evenly distributed from the raw data. For each sampling point, the processor extracts its X-axis, Y-axis, and Z-axis coordinates in the camera coordinate system, as well as the corresponding polarized blue light grayscale value after distortion correction.

[0028] Furthermore, to enhance the model's ability to recognize surface geometry, the processor calculates a normal vector (Nx, Ny, Nz) based on the 20 nearest neighbors of each sampling point. The processor concatenates the coordinate values ​​(3D), the normal vector (3D), and the polarization grayscale values ​​(1D) to generate a feature tensor with a dimension of 4096×7. To eliminate dimensional differences, the coordinate values ​​are normalized to the interval [-1, 1], and the polarization grayscale values ​​are normalized to the interval [0, 1]. This feature tensor is used as input data and transmitted to the subsequent semantic segmentation neural network model.

[0029] The semantic segmentation neural network model adopts an architecture that combines a point-based multilayer perceptron with local feature aggregation: Feature encoding layer: 1×1 convolution, output channels 64 and 128, activation function ReLU; local topological features: K-nearest neighbor search (k=20), calculate relative coordinates and concatenate with neighborhood features; surface attribute fusion: polarization grayscale is fused with geometric features through a weight matrix; global feature aggregation: max pooling generates a 1024-dimensional global feature vector; decoding layer: fully connected layer (512, 256, 128 nodes), dropout rate 0.5; output layer: 4096×3 matrix, representing three state categories (rigid, flexible, background). Classification and Mask Generation: Softmax normalization is used to calculate the probability distribution, and state labels are assigned based on the maximum posterior probability; a binary Boolean classification mask (True / False) of length 4096 is generated, and the corresponding subspace is extracted using logical indexes.

[0030] Furthermore, regarding the architecture and feature extraction of the semantic segmentation neural network model, specifically, the semantic segmentation neural network model adopts an architecture combining a point-based multilayer perceptron (MLP) with a local feature aggregation module; the specific execution flow of this model is as follows: Local topological feature extraction: The input tensor first enters the first-layer feature encoding module; this module contains two consecutive convolutional layers with a kernel size of 1×1 (i.e., processing each point independently), and the number of output channels are set to 64 and 128 respectively, with ReLU activation function; in order to capture local topological features, before each convolutional layer, the model performs a K-nearest neighbor (k-NN) search, with k=20; for each center point, the system finds its 20 nearest Euclidean neighbors, calculates the relative coordinate difference between the neighbors and the center point, and concatenates the relative coordinate difference with the feature vector of the neighbors themselves; this processing method can explicitly encode the geometric microstructure of local regions, such as edges, corners, or curvature changes, thereby extracting local topological features.

[0031] Surface property feature fusion: Since the input data explicitly contains polarized blue light grayscale values ​​(surface optical property information), this information participates in the computation as an independent channel in the network; in the middle layer of the network (i.e., the feature transformation layer that maps from 128 dimensions to 256 dimensions), the polarization grayscale features are weighted and fused with the geometric features through the multiplication operation of the weight matrix; since the polarization reflection characteristics of metal surfaces are very different from those of rubber or plastic, this step enables the network to learn the nonlinear relationship between material properties and geometry, and extract high-dimensional surface property features.

[0032] Global Feature Aggregation and Classification Probability Calculation: After multi-layer feature extraction, the data dimension becomes 4096×1024. At this point, the processor performs a max pooling operation, selecting the feature with the largest response value in the feature channel dimension to generate a 1024-dimensional global feature vector. Subsequently, this global feature vector is copied and concatenated with the local features of each point, and decoded through three fully connected layers (with 512, 256, and 128 nodes respectively). Each layer is configured with Dropout (dropout rate set to 0.5) to prevent overfitting. Finally, the output layer is a matrix with a dimension of 4096×3, where "3" represents three preset state categories: rigid state (such as metal parts), flexible state (such as cables / hose), and background noise.

[0033] Furthermore, the calculation of state class probabilities and the determination of maximum a posteriori probabilities are specifically addressed by the fact that the process of calculating the state class probability for each data point is not a black-box output, but rather implemented through a standardized exponential function (Softmax); for the i-th data point, the three raw scores are output. The processor calculates its natural index respectively. Each index value is then divided by the sum of the three index values. For example, the logic for calculating the probability of a rigid state is: the index value of the rigid fraction is divided by (the index value of the rigid fraction + the index value of the flexible fraction + the index value of the background fraction). This step ensures that the sum of the output probability values ​​is 1 and that the value is between 0 and 1.

[0034] Subsequently, the processor performs classification using the Maximum A posteriori (MAP) principle. The system's decision logic is as follows: For each data point, its probability values ​​in the three dimensions of "rigidity," "flexibility," and "background" are compared. If the rigidity probability of a point is 0.1, the flexibility probability is 0.85, and the background probability is 0.05, the processor determines that the posterior probability of that point is the highest in the flexible category, and thus assigns the "flexible" state label to that data point. This process iterates through all 4096 points, generating a state label vector of length 4096.

[0035] Furthermore, the classification mask generation and spatial filtering specifically involve the processor generating rigid and flexible classification masks based on the generated state label vectors; the masks are essentially Boolean index arrays. The rigid classification mask is a binary array of length 4096. The value of the corresponding position is True (or 1) if and only if the label of a point is "rigid", otherwise it is False (or 0). The flexible classification mask is similar. The value of the corresponding position is True if and only if the label of a point is "flexible".

[0036] Finally, the classification mask is applied to perform spatial filtering on the original observation data; the processor uses logical indexing technology to perform dot product operation (or index extraction) on the original 4096×7 data tensor and the rigid classification mask; all rows with a mask value of True are retained and recombined to construct a rigid state subspace; similarly, a flexible state subspace is extracted using a flexible classification mask.

[0037] By introducing a semantic segmentation neural network to construct a rigid-flexible classification mask, the observation data is accurately isolated at the feature level, avoiding interference from background noise or heterogeneous objects on the solution of specific states. This spatial filtering mechanism significantly improves the accuracy of target recognition in complex stacking scenarios.

[0038] Furthermore, the calculation of feature alignment and positional deviation in the rigid state subspace is addressed. Specifically, the processing of the rigid state subspace aims to solve the problem of high-precision positioning of the part pose; the input of this step is the filtered rigid point cloud data (source point set), and the preset data is the rigid body standard model point cloud (target point set, usually generated by discretization of CAD model, with approximately 2000 points) stored in the database; the processor uses an improved Iterative Closest Point (ICP) registration algorithm to perform feature alignment. This process does not directly call black-box functions, but includes the following specific closed-loop calculation logic: Spatial Index Construction and Corresponding Point Search: To accelerate the search process, the processor first constructs a KD-tree spatial index structure for the target point set. At the beginning of each iteration, for each 3D coordinate point in the source point set, the processor uses the KD-tree to search for its nearest Euclidean neighbor in the target point set. To eliminate noise or interference from non-overlapping regions, the system sets a distance threshold of 5.0 mm. If the calculated distance between point pairs is greater than this threshold, the point pair is determined to be an invalid match and will not participate in the subsequent transformation matrix calculation.

[0039] The transformation matrix is ​​solved (based on SVD decomposition): Based on the selected set of valid point pairs, the processor first calculates the centroid of the source point set and the centroid of the target point set respectively; then, a decentralization operation is performed, that is, the coordinate values ​​of all points are subtracted from their corresponding centroid coordinates; the processor uses the transpose of the decentralized source point set matrix and the target point set matrix to perform matrix multiplication to construct the covariance matrix (3×3 dimension).

[0040] Furthermore, singular value decomposition (SVD) is performed on the covariance matrix to decompose it into a left singular vector matrix U and a right singular vector matrix V; the rotation matrix for the current iteration step is obtained by calculating the product of matrix V and the transpose of matrix U; the translation vector is obtained by subtracting the source centroid after rotation transformation from the target centroid. The system's default convergence conditions are: the change in root mean square error (RMSE) is less than 0.01 mm, or the maximum number of iterations reaches 50. After each calculation of the transformation parameters, the processor applies them to the source point set and calculates the new average distance between point pairs. If the convergence conditions are not met, the system uses the updated point set to enter the next iteration. If the conditions are met, the system outputs the final cumulative transformation matrix (including the final rotation matrix and translation vector).

[0041] Position deviation matrix generation: The processor retrieves the pre-set ideal working coordinate system matrix (representing the standard theoretical pose of the robot end effector when grasping the part); the calculation logic of the position deviation matrix is: using the inverse matrix of the ideal working coordinate system matrix, multiplying it by the final cumulative transformation matrix output by the ICP algorithm; this deviation matrix accurately describes the six degrees of freedom offset of the current part relative to the standard grasping pose (X, Y, Z axis displacement and roll, pitch, and yaw angle deviations), which serves as the compensation input for subsequent robot path planning.

[0042] Furthermore, geometric topological reconstruction and deformation analysis of flexible state subspaces are conducted. Specifically, for the flexible state subspace (typically cable data distributed in tubular or strip shapes), rigid registration cannot be used due to its non-fixed shape. The processor employs an elliptical cylindrical manifold fitting algorithm to analyze its morphology. The specific steps are as follows: Principal axis direction vector extraction (PCA analysis): The processor first performs principal component analysis (PCA) on all point cloud data within the flexible subspace. Specifically, the processor calculates the covariance matrix of the point cloud coordinate data and performs eigenvalue decomposition on the matrix. The system selects the eigenvector corresponding to the largest eigenvalue as the principal axis direction vector of the flexible operation object. This vector physically represents the overall spatial extension trend of the cable.

[0043] Cross-sectional projection and topology reconstruction: The processor constructs a two-dimensional projection plane with a normal vector parallel to the aforementioned principal axis direction vector; all three-dimensional points in the flexible subspace are projected onto this plane along the principal axis direction to form a two-dimensional cross-sectional scatter plot; subsequently, the least squares method based on Random Sample Consensus (RANSAC) is used to fit the two-dimensional ellipse equation.

[0044] The specific calculation logic is as follows: The system randomly selects 5 data points from the projection point set and solves the coefficients of the ellipse equation by minimizing the algebraic distance; then it calculates the geometric distance from all other points to the fitted ellipse. If the distance is less than 0.5 mm, the point is marked as an interior point. This process is repeated 100 times, and finally the model parameter with the highest proportion of interior points is selected as the optimal fitting result, thereby reconstructing the surface geometric topology of the flexible object.

[0045] Geometric feature parameter calculation: Based on the best-fit ellipse geometric parameters, the processor extracts the lengths of the major and minor semi-axis. Major-minor axis ratio calculation: The processor divides the value of the major semi-axis by the value of the minor semi-axis; this ratio reflects the degree of flatness of the cable cross-section (for example, the ratio of a round cable should be close to 1.0). Spindle tilt angle parameter calculation: The processor calculates the cosine value of the angle between the spindle direction vector and the Z-axis (vertical axis) of the robot base coordinate system, and obtains the specific tilt angle value (in degrees) through the inverse cosine function. Deformation state variable generation: The system compares the calculated real-time parameters with the preset standard object shape data; the preset data includes the standard cross-sectional ratio (e.g., set to 1.05, allowing for minor tolerances) and the standard installation tilt angle (e.g., set to 0 degrees, i.e., vertical installation). The processor calculates the absolute value of the difference between the major and minor axis ratios and the absolute value of the difference between the principal axis tilt angles through subtraction operations; these two differences are combined into a two-dimensional vector, namely the deformation state variable.

[0046] Furthermore, to verify the feasibility of the above calculations, the processor sets a boundary check before performing manifold fitting. If the number of valid points in the flexible subspace is less than 50 (e.g., due to data loss caused by occlusion), the processor will not perform the fitting operation, but will output the deformation state variable as an invalid value (e.g., -1) and trigger an anomaly signal. Finally, the calculated position deviation matrix (from the rigid processing path) and deformation state variables (from the flexible processing path) are packaged into a unified state data packet. This data packet is sent to the robot motion controller in real time via the industrial bus at a period of 10ms. The controller corrects the end coordinates of the robotic arm based on the position deviation matrix, and at the same time determines whether the gripping force needs to be adjusted based on the deformation state variables (for example, if the ratio of major axis to minor axis is abnormally large, it means that the cable is squeezed and the gripper pressure needs to be reduced), thereby realizing adaptive operation based on physical state perception.

[0047] Through the above steps, the system successfully decomposes the mixed observation data into two independent subsets of physical properties at the data stream level; the rigid subspace data will be sent to the subsequent pose estimation algorithm (such as ICP registration), while the flexible subspace data will be sent to a dedicated deformation analysis module, thereby realizing differentiated processing for heterogeneous objects.

[0048] Furthermore, the capture and construction of the time-domain response sequence, specifically, the step of capturing the horizontal vibration displacement of the top of the column and constructing the time-domain response sequence, is completed by a high-precision laser displacement sensor deployed at the end of the actuator in conjunction with a high-speed data acquisition card.

[0049] During the pre-calibration phase before the start of the contact operation, the system applies a standard pulse excitation force to the column (e.g., an instantaneous thrust with a duration of 50ms and an amplitude of 10N applied by the end of the robotic arm), inducing the column to generate free decaying vibration; at this time, the data acquisition card starts acquisition at a sampling frequency of 1000Hz, which is set based on the Nyquist sampling theorem and is sufficient to cover the estimated low-frequency (usually below 100Hz) structural vibration characteristics of the column. The acquisition process lasted for two seconds, acquiring two thousand discrete displacement sampling points. The processor linearized the original voltage signal to obtain physical displacement values ​​in millimeters. To eliminate interference from high-frequency electromagnetic noise in the environment, the processor performed a moving average filter on the original data, setting the sliding window size to 5, i.e., taking the average of the current point and the two points before and after it as valid data, ultimately generating a smooth "time-displacement" time-domain response sequence, which served as input data for subsequent spectrum analysis. Differentiated solution strategies were adopted for subspaces with different physical characteristics, ensuring both the convergence accuracy of the pose calculation of rigid components and the precise quantification of the geometric deformation of flexible cargo. This provided accurate position deviation feedback and deformation state indicators for the subsequent control system, achieving precise perception with "targeted treatment".

[0050] Furthermore, for the above time-domain response sequence, the processor performs a Fast Fourier Transform (FFT) to extract the first-order natural frequency; In terms of computational logic, the processor first performs zero-padding on the 2000-point time-domain sequence, extending its length to 2048 points (i.e., ...). To meet the computational requirements of the radix-2 FFT algorithm, the processor applies Hanning window weighting to the time-domain data to reduce spectral leakage. After the FFT operation, the system obtains a frequency domain array containing amplitude and phase information. The processor traverses the amplitude spectrum of this array and searches for the spectral peak with the largest amplitude within the effective frequency band of 0.5-50Hz. The frequency index value corresponding to this spectral peak is the first-order natural frequency of the pillar. For example, in a certain measurement, the calculated main peak frequency was 12.5Hz. Subsequently, the system uses the logarithmic decay method to estimate the damping ratio; the processor returns to the time-domain response sequence and uses a peak search algorithm to identify two adjacent positive maxima in the vibration waveform (denoted as the first peak and the second peak); the calculation logic is as follows: take the ratio of the first peak to the second peak, calculate the natural logarithm of the ratio, and then divide it by twice pi (approximately 6.28) to obtain the equivalent damping ratio of the system; for example, if the first peak is measured to be 1.5 mm and the second peak is 1.2 mm, the damping ratio can be calculated to be approximately 0.035; this parameter quantifies the energy dissipation capacity of the column itself and provides a physical benchmark for tuning the damping coefficient in the subsequent impedance model.

[0051] Furthermore, the establishment and parameter tuning of the impedance contact model, specifically, the establishment of a second-order linear differential equation describing the dynamic mapping between the position correction and the contact force as the impedance contact model, is the core of achieving compliant control. The model is discretized into a difference equation within the controller for numerical computation; the model accepts contact force deviation (i.e., the difference between the measured value of the force sensor and the expected contact force) as input and outputs the position correction of the end effector; the model contains three key physical parameters: virtual mass, damping coefficient, and stiffness coefficient.

[0052] Virtual mass setting: characterizes the system's inertial response to sudden changes in contact force; the system sets it to 5KG; the larger the value, the slower and smoother the system's response to force impact; the smaller the value, the more sensitive the response, but it may cause high-frequency oscillations. Damping coefficient setting: characterizes the rate at which the system dissipates vibration energy; the processor calculates the reference damping value based on the actual damping ratio (0.035) and first natural frequency (12.5Hz) of the column identified in the previous steps, combined with the critical damping formula, and introduces a safety factor of 1.2 times, setting the damping coefficient to 100 N / m; this ensures that the robotic arm can quickly suppress vibration during contact and avoid continuous contact jitter.

[0053] Stiffness coefficient setting: characterizes the elastic restoring force of the system to recover to the desired trajectory after being subjected to external force; the system sets it to 2000 N / m; By using frequency domain analysis and logarithmic decay method to obtain the time-varying first-order natural frequency and damping ratio in real time, the model parameter mismatch problem caused by height increase and load change of stacker crane is solved. The established second-order impedance model accurately maps the dynamic contact characteristics of the system, laying the physical model foundation for achieving stable compliant control.

[0054] At each time step of the control cycle (e.g., every 1ms), the processor executes the following computational logic: First, the current contact force deviation is calculated. Second, using the position correction, velocity correction, and current force deviation from the previous moment, and based on the virtual mass, damping coefficient, and stiffness coefficient set above, the second-order differential equation is solved using numerical integration methods (such as trapezoidal integration or Runge-Kutta method). Finally, the position correction for the current moment is output. For example, when an unexpected impact force of five Newtons is detected, the impedance model will calculate a small backward displacement (such as 0.05mm), which is executed through the position loop command of the robotic arm controller, thereby flexibly absorbing the impact energy, maintaining the stability of the contact force, and realizing dynamic closed-loop adjustment from "force input" to "displacement output".

[0055] Furthermore, based on the construction of the desired velocity vector for visual servoing, specifically, the impedance control module first receives the position deviation matrix (including displacement deviations of the X, Y, and Z axes and angular deviations of the three attitude angles) from the geometric analysis module; in order to convert these deviations in Cartesian space into motion commands for each joint of the robotic arm, the processor adopts a position-based visual servoing control law.

[0056] First, the processor constructs the Jacobian matrix; this matrix is ​​a real-time transformation matrix with a dimension of 6×N (N is the number of axes of the robotic arm, for example, 6 axes), and its value is obtained by calculating the angle values ​​fed back by the encoders of each joint at the current moment through the positive kinematics differential, which represents the linear relationship between the joint space velocity and the Cartesian space velocity. Next, the processor performs inverse kinematics calculation and proportional gain adjustment. The calculation logic is not a simple matrix inversion, but uses the damped least squares method to solve the pseudo-inverse of the Jacobian matrix to avoid numerical instability under singular configurations. The processor takes the six components in the position deviation matrix as error inputs and multiplies them by preset proportional gain coefficients. For example, the position gain is set to 2.5 and the attitude gain to 3.0. This gain value determines the robot's response speed to eliminate errors. The calculated Cartesian space expected velocity (error multiplied by gain) is then left-multiplied by the Jacobian pseudo-inverse matrix to solve for the expected velocity vector of each servo axis, i.e., the angular velocity command of the six joints. This vector serves as the system's reference input and does not yet include compliance characteristics.

[0057] Furthermore, to endow the robot with adaptability to flexible objects, the system dynamically adjusts the parameters of the impedance model based on the deformation state variables output from previous steps (e.g., the deviation of the major and minor axis ratio of the cable cross-section). The deformation state variable Δs = [Δr, Δθ]^T is a two-dimensional vector, where Δr is the deviation of the major and minor axis ratio from the standard value (dimensionless), and Δθ is the deviation of the principal axis tilt angle from the standard tilt angle (in degrees or radians). The processor establishes a nonlinear attenuation mapping function between the deformation state variables and the stiffness coefficient. This function is physically designed as a "softening protection mechanism": when an increase in object deformation is detected, the stiffness of the robotic arm's end effector is actively reduced to decrease the contact force. The specific algorithm logic is as follows: The basic stiffness coefficient is set at 2000 N / m, corresponding to the rigid working state without deformation, and the minimum stiffness limit is set at 500 N / m. The mapping function adopts the inverse S-shaped (Sigmoid) curve logic. For example, when the deformation state variable (ratio deviation) is in the range of 0 to 0.1, the stiffness remains unchanged at 2000 N / m. When the deviation exceeds 0.1 and increases towards 0.5, the stiffness coefficient decreases rapidly according to the inverse square relationship of the deviation value. When the deviation reaches 0.5, the stiffness coefficient decreases to 500 N / m.

[0058] Based on the constant damping ratio constraint, the system synchronously adjusts the damping coefficient. To ensure that the system remains in a critically damped or overdamped state during stiffness changes, and to avoid position overshoot or oscillation, the system sets a constant target damping ratio of 0.85. Within each control cycle, e.g., 1 ms, the processor first obtains the current real-time stiffness coefficient and a preset virtual mass, e.g., 5 kg, and then sets the constant target damping ratio. For example, 0.85; the processor adjusts the stiffness coefficient according to the real-time changing stiffness coefficient. and preset virtual quality Calculate the damping system that should be applied at the current moment. For example, when the stiffness decreases from 2000 to 500, the damping coefficient will automatically decrease from about 170 N / m to about 85 N / m to ensure the consistency of dynamic response characteristics. Furthermore, the execution and command generation of the admittance control structure, specifically the process of generating the final control command using the admittance control structure, involves "dynamic filtering" of the desired velocity vector. The system uses the desired velocity vector calculated in the previous step as a reference input, and substitutes the real-time adjusted stiffness coefficient, such as 800 N / m, damping coefficient, such as 110 N / m, and virtual mass of 5 kg into the second-order impedance contact model to construct a digital dynamic filter. The execution logic of the dynamic filter uses the second-order Runge-Kutta method for discrete-time integration; within each 1ms control cycle, the processor performs the following operations: The driving force is calculated by multiplying the difference between the desired velocity vector and the actual feedback velocity vector by the damping coefficient, and then adding the difference between the desired position and the actual position by the stiffness coefficient to obtain the virtual "corrective force". The acceleration is obtained by dividing the above "corrective force" by the virtual mass. The integral solution involves the processor using Euler's integral method to convert acceleration into velocity change, updating and correcting the velocity, and then performing the integral operation again to convert it into position change, thus obtaining the position correction. Finally, the processor superimposes the position correction onto the original trajectory planning instruction to generate the initial control instruction. This instruction is then sent to the servo drive via the EtherCAT industrial bus.

[0059] Example of a closed-loop data flow: Suppose the visual servo calculates that the robotic arm should move forward at a speed of 10 millimeters per second (desired velocity vector); at this time, the flexible sensing module detects that the cable is slightly compressed, with a deformation variable of 0.2. The parameter adjustment module reduces the stiffness coefficient from 2000 to 1200 N / m through a mapping function, and the damping coefficient is adjusted synchronously. The admittance filter receives an input of 10 mm / s, but due to the reduction in stiffness, the model exhibits greater compliance. If contact resistance is encountered at this time, the position correction calculated by the filter will slow the actual forward movement speed of the robotic arm to 8 mm / s and generate a 0.5 mm backward compensation in position. The final control command does not force a 10 mm / s advance, but includes a 0.5 mm flexible avoidance, thus ensuring proper assembly while physically protecting the flexible cable from being crushed.

[0060] A deformation-based variable stiffness admittance control strategy was constructed, which can automatically reduce the system stiffness to achieve compliant contact when the cargo deforms, effectively preventing cargo damage. Simultaneously, combined with a visual servo closed-loop, the end effector achieves precise tracking and correction of the desired trajectory while ensuring flexible interaction.

[0061] Furthermore, a multi-rate Kalman filter model and prior state estimation are constructed. Specifically, the estimation and shaping module first establishes a kinematic recursive model for processing heterogeneous sensor data; since the sampling frequency of the robot joint servo encoder is much higher than that of the vision sensor, the system adopts a multi-rate Kalman filter architecture; Data source and sampling settings: The system sets the sampling frequency of the servo encoder to 1000Hz as a high-frequency input to provide joint angle and angular velocity data. The system sets the sampling frequency of the vision sensor to 30Hz as a low-frequency observation input, providing a position deviation matrix in Cartesian space; Kinematic recursion and time update: Within each 1ms control cycle, the processor first uses the forward kinematics equations to convert the joint angles fed back by the encoder into the position vector of the end effector in the Cartesian coordinate system. and velocity vector The processor constructs a state vector X, which contains position and velocity components corresponding to the six degrees of freedom (dimension 12×1). The time update (prediction) steps for performing Kalman filtering are as follows: The processor uses a constant velocity motion model as the state transition matrix A (where the predicted position is equal to the position at the previous moment plus the velocity at the previous moment multiplied by the time step of 0.001 seconds). Calculate the prior state estimate That is, the current state calculated based on the optimal estimate of k from the previous time step; At the same time, the processor updates the covariance matrix. The calculation involves multiplication and addition operations of the state transition matrix A and the preset process noise matrix Q (the diagonal elements are set to 1e-5 to represent the uncertainty of the model); Furthermore, when the timeline reaches the visual sampling point, the system performs a measurement update step to eliminate accumulated drift; Observation vector construction and residual calculation: The processor receives the position deviation matrix from the vision module and converts it into observation vectors. ; Computational visual observation residuals: ;in This represents the residual vector at time K. H represents the visual observation vector at the current moment, and H is the measurement matrix used to extract the position vector from the state vector. This represents the prior state estimate.

[0062] Furthermore, the processor calculates the Kalman gain, a parameter that determines the level of trust in the visual observation data; the system presets the observation noise covariance matrix R, with its diagonal elements set to 1e-3, representing the noise level of the visual measurement; The calculation formula is as follows: ;in This represents the optimal Kalman gain matrix. R represents the prediction error covariance matrix, and R represents the system preset observation noise covariance matrix. This indicates the transpose of the measurement matrix; The processor uses the calculated gain to correct the system state, incorporates visual feedback information into the state variables, eliminates cumulative drift, and obtains the optimal posterior state estimate; the calculation formula is as follows: ;in, This represents the corrected posterior state estimate (including optimal position and velocity). This represents the prior state estimate before correction; Finally, based on the corrected optimal position, the processor generates speed control commands to drive the servo motor; the system sets the position loop proportional gain. The value is set to 50.0 to control the response speed for error elimination. The specific calculation formula is as follows: ,in This indicates the generated corrected speed control command. For position loop proportional gain, This represents the theoretical target position output by the trajectory planning module. Represents the estimated value of the posterior state. The optimal position component extracted from it; Furthermore, the construction of the zero-vibration pulse sequence, specifically, in order to eliminate the residual vibration generated when the robotic arm stops moving at high speed, the processor constructs an input shaper based on the physical parameters identified by the preceding dynamic modeling module; Parameter extraction: The processor reads the first-order natural frequency fn (e.g., 12.5) and damping ratio ζ (e.g., 0.035) from memory. Pulse time interval calculation: The system adopts a zero-vibration (ZV) shaping strategy, which decomposes a complete input instruction into two time-staggered pulses; Calculation of vibration damping period Td: The processor first calculates the undamped natural period (0.08 seconds) using the formula 1 / fn, and then fine-tunes it according to the damping ratio (since the damping ratio is very small, the damping period is approximately 0.08 seconds). Calculate the pulse time interval Δt: set it to half of the damping period, i.e., Δt = Td / 2; combining the above values, we calculate Δt = 40ms; this means that the second pulse will be emitted 40ms later than the first pulse to produce destructive interference; Amplitude coefficient calculation: The processor calculates the pulse amplitude coefficient using the damping ratio ζ; First, calculate the attenuation factor. Substituting ζ=0.035, we calculate K≈0.895; calculate the amplitude of the first pulse. Calculate the amplitude of the second pulse. The zero-vibration pulse sequence constructed thus is as follows: the amplitude is 0.528 at t=0 and the amplitude is 0.472 at t=40ms. Furthermore, the temporal convolution operation and the generation of shaping execution instructions are specifically: the processor performs convolution operations in the temporal domain to transform the aforementioned modified speed control instructions into the final execution instructions; Instruction buffer queue management: The processor allocates a 50-byte first-in-first-out (FIFO) circular buffer in memory (corresponding to 50ms of historical data, which is sufficient to cover the 40ms latency requirement); the "corrected speed control instruction" generated at each moment is pushed into this buffer.

[0063] Real-time convolution calculation: Within each 1ms interpolation cycle, the processor performs the following weighted summation operation: take the original instruction value at the current time t and multiply it by the first pulse amplitude A1 (0.528); access the buffer, read the historical instruction value at time t−40ms and multiply it by the second pulse amplitude A2 (0.472); add the products of the above two parts to obtain the shaped execution instruction at the current time; after convolution processing, the original step speed instruction is transformed into a two-level stepped pulse sequence; physically, the first component drives the robotic arm to move and excites vibration, and the second component, which arrives 40ms later, drives the robotic arm to continue moving and excites a vibration wave with the same amplitude and opposite phase as the previous vibration, thereby physically canceling the residual sway of the column; Ultimately, the shaping execution command is converted into a PWM signal or bus command and sent to the servo driver to drive the motor to perform vibration-free high-speed operation.

[0064] By fusing vision and encoder data using multi-rate Kalman filtering, the limitations of single-sensor bandwidth and noise issues are resolved, resulting in high-frequency and smooth control commands. Combined with zero-vibration input shaping technology, the identified dynamic parameters are used to cancel residual oscillations in the mechanical structure at the command source, achieving rapid, overshoot-free positioning of the system.

[0065] Furthermore, the construction and boundary definition of the virtual repulsive potential field are specifically achieved by the processor first receiving surface geometric topology data output by the preceding geometric analysis module; in stacker crane applications, this data is typically represented as fitted cargo outline parameters, such as an elliptical cylindrical equation describing the bulging surface of a cardboard box. The adaptive potential field boundary control method based on elliptical cylindrical manifolds first performs surface geometric topology analysis on the target object, and then establishes its mathematical model by fitting an elliptical cylindrical manifold. The system calculates key geometric feature parameters, specifically including the lengths of the major and minor semi-axes of the elliptical cross-section, the spatial center coordinates of the cylinder, and the direction vector describing the orientation of the cylinder's axis. Based on these features, a parametric model of the elliptical cylinder, determined by both the circumferential angle and the axial height, is constructed.

[0066] To improve the efficiency of real-time computation, the aforementioned parametric model is discretized. The number of sampling points is set in both the circumferential and axial dimensions (e.g., 36 points circumferentially and 8 points axially), generating hundreds of mesh vertices covering the surface. The spatial data of these vertices is stored in a lookup table for subsequent fast indexing and matching.

[0067] During system operation, the current spatial position of the end effector is acquired in real time. The algorithm traverses a lookup table to calculate the Euclidean distance from the end effector position to each grid vertex. Through traversal and comparison, the algorithm selects the minimum distance value and the nearest point on the surface corresponding to that distance. To determine the correct direction of the obstacle avoidance repulsion force, the algorithm analyzes the selected nearest point and its surrounding neighborhood points. First, the tangent vectors of the surface in two orthogonal directions are calculated, and then the normal vector perpendicular to the surface is obtained through vector cross product. Finally, this normal vector is normalized to obtain the unit normal vector indicating the outward direction of the surface.

[0068] The system pre-sets a safe influence range (e.g., 50 mm) and introduces deformation state variables and corresponding gain coefficients to achieve adaptive adjustment.

[0069] When the minimum distance between the end effector and the surface is detected to be less than the set safety influence range, the virtual repulsion mechanism is triggered: Amplitude calculation is performed based on the depth of the intrusion into the safe zone to determine the magnitude of the repulsion velocity. The closer the distance, the stronger the repulsion, and this strength is adjusted by the deformation state gain. Direction determination involves setting the direction of the repulsion velocity to the opposite direction of the unit normal vector of the surface. The calculated virtual repulsion velocity is vector-superimposed with the system's desired reference velocity to output the final composite control velocity, thereby achieving compliant avoidance of obstacle boundaries while ensuring the execution of the mission path.

[0070] The processor constructs a virtual repulsive potential field based on these geometric parameters; the system defines the maximum contour envelope of the flexible cargo (i.e., the geometric surface that protrudes outward most from the cargo surface) as the zero-distance boundary of the potential field; in order to achieve a balance between operational efficiency and safety, the system sets a potential field influence range at the software level; for example, this range is set as the area extending outward 50mm from the zero-distance boundary in the normal direction.

[0071] In this step, the data flow logic is as follows: the input is the major semi-axis, minor semi-axis and center coordinates of the elliptical cylinder; the output is the three-dimensional spatial domain of the potential field; within this domain, each coordinate point in the space is assigned a non-zero potential field strength value, while outside the domain (i.e., more than 50mm away from the goods), the potential field strength is 0 and has no effect on the stacker crane.

[0072] Furthermore, in order to achieve adaptive control, the processor establishes a positive correlation between the potential field intensity gradient and the deformation state variables (such as the major and minor axis ratio deviation values ​​calculated in the previous step). The purpose of this step is that when the cargo is more severely deformed (indicating that the material is more fragile or the internal pressure is greater), the "repulsive force" generated by the virtual potential field should be stronger, thereby forcing the stacker crane to decelerate earlier and more gently. The specific algorithm logic adopts piecewise linear gain adjustment: the processor sets the base gain coefficient to 1.0; the threshold of the deformation state variable is set to 0.2; within each 10ms control cycle, the processor reads the current deformation state variable; if the variable is less than or equal to 0.2, the potential field strength gain remains at the base value (1.0); if the variable is greater than 0.2, the processor performs linear amplification calculation: subtract the threshold from the current variable value, multiply by the preset sensitivity coefficient (e.g., set to 10), and then add the result to the base value; for example, when the deformation variable reaches 0.3, the calculated potential field strength gain is 2.0; this means that at the same distance, the system will produce twice the avoidance response of the normal situation.

[0073] Furthermore, Euclidean distance gradient calculation and vector generation, specifically, during the generation of initial control commands, the processor monitors in real time the positional relationship between the stacker crane end effector (such as the fork tip) and the potential field boundary; Distance calculation: The processor obtains the Cartesian space coordinates of the end effector at the current moment (obtained by forward kinematics calculation from encoder feedback); the processor calculates the shortest Euclidean distance from this coordinate point to the aforementioned "zero-distance boundary" (i.e., the elliptical cylinder); at the same time, the processor calculates the unit normal vector from the boundary surface to the end effector, which represents the direction of the repulsive force. Virtual repulsion velocity vector generation: If the calculated Euclidean distance is less than the previously set potential field influence range (50mm), the algorithm triggers the repulsion logic; The processor uses an inverse proportional function to generate the repulsive velocity amplitude: the difference between (potential field influence range minus the current Euclidean distance) is divided by the potential field influence range to obtain a normalization factor between 0 and 1; Then, the normalization factor is multiplied by the maximum allowable repulsion velocity (e.g., set to 100 mm / s), and then multiplied by the potential field strength gain calculated in step two (e.g., 2.0). Finally, the calculated amplitude is assigned to the aforementioned unit normal vector to generate a virtual repulsive velocity vector in three-dimensional space; For example, when the fork tip is only 10mm away from the surface of the goods, and the goods are deformed to a large extent resulting in a gain of 2.0, the system calculation logic is as follows: normalization factor = (50-10) / 50 = 0.8; at this time, the repulsion velocity amplitude = 0.8 × 100mm / s × 2.0 = 160mm / s; the system generates a strong repulsion velocity of 160mm / s along the normal direction backward. Furthermore, the speed superposition and reverse damping effect are achieved. Specifically, the processor performs a superposition operation on the speed field to generate the final input instruction to the servo system. The processor acquires the desired velocity vector generated by the trajectory planning module (e.g., the forks are extending forward at a speed of 200 mm / s).

[0074] The processor performs a vector addition operation between the virtual repulsion velocity vector calculated above and the desired velocity vector; Since the direction of the virtual repulsion velocity is usually opposite to the direction of the fork's movement toward the cargo, this vector superposition physically manifests as a virtual reverse damping force. Continuing with the above data: if the desired forward velocity is 200 mm / s, and the virtual repulsion velocity is 160 mm / s backward, the net velocity after superposition is only 40 mm / s; This not only significantly reduces the approach speed, but also, if the forks continue to approach, the repulsion speed will increase rapidly and eventually exceed the desired speed, causing the net speed to become negative (i.e., automatically retract).

[0075] Finally, the superimposed composite velocity vector is fed into the impedance control model (as a reference input), and after being processed by the admittance filter, it is converted into torque or speed commands for each joint motor. This process ensures that the stacker crane can "sense" the danger based on the degree of deformation of the goods before it comes into contact with them, and automatically generate a flexible deceleration and hovering effect similar to the repulsion of like poles of magnets, thereby completely avoiding the damage of rigid impacts to fragile goods. A virtual potential field guidance mechanism based on deformation boundaries is introduced, adding an active safety barrier to the operation process. This mechanism generates a reverse repulsion velocity during extreme deformation or approach to the boundary. Combined with impedance control, this further mitigates the risk of excessive compression, significantly improving the system's obstacle avoidance capability and safety in complex or confined space operations.

[0076] This embodiment achieves high-precision pose control and adaptive compliant interaction for large-inertia stacker cranes by decoupling observation of rigid-flexible states and online identification of dynamic parameters. Utilizing the combined effects of impedance parameter modulation and input shaping techniques, residual vibration of the columns during variable-speed motion is effectively eliminated, significantly improving the efficiency of automated warehousing operations and cargo safety. Example

[0077] This embodiment focuses on the dynamic stacking operation of liquid chemical soft bags stored in a large automated warehouse. Such goods not only have flexible and easily deformable surfaces, but also the internal liquid sloshing will generate nonlinear time-varying interference forces, which puts forward higher requirements for the predictive ability and response speed of the control system. The embodiment focuses on describing two key steps: "deformation trend prediction based on time-series deep learning" and "impedance parameter self-tuning based on multi-objective optimization". First, based on the deformation trend prediction of the LSTM network, specifically, in order to cope with the hysteresis deformation effect of liquid soft package during the transportation process, the system does not rely solely on the current instantaneous observation value, but introduces a Long Short-Term Memory (LSTM) network to predict the future deformation state in advance. Further, the construction of the temporal feature tensor involves the processor first constructing a time sliding window with a length of 50 control cycles (covering historical data over the past 500ms). For each time step within the sliding window, the system collects physical quantities in four key dimensions: the ratio of the major and minor axes of the flexible cargo, the tilt angle of the main axis, the vertical loading speed of the end effector, and the feedback value of the contact force sensor. The processor standardizes these 50 sets of data to construct a temporal feature tensor with a dimension of 50×4. The construction and preprocessing of the feature tensor involve the processor performing farthest-point sampling on the original observation data, selecting 4096 representative data points. For each sampling point, coordinate values ​​(3D), normal vectors (3D), and polarization grayscale values ​​(1D) are extracted to construct a feature tensor with a dimension of 4096×7. The coordinate components are normalized to the [-1,1] interval, and the grayscale components to the [0,1] interval. This serves as the input to the deep learning model. The architecture and execution of the prediction model are as follows: This input tensor is fed into a double-layer stacked LSTM neural network model. The first layer is the sequence processing layer, which contains 64 hidden layer nodes. The activation function is the hyperbolic tangent function (Tanh), which is used to extract long-term dependency features in the time series and output sequence vectors of the same dimension. The second layer is the state compression layer, which also contains 64 hidden layer nodes, used to fuse features and output the hidden state vector of the last time step; The model's output is connected to a fully connected layer containing one neuron, using a linear activation function; In each inference process, the model outputs the estimated deformation state variable value for the next 10ms based on the current and past 500ms data. If the model predicts that the future deformation variable value exceeds the preset safety threshold (e.g., 0.3), the system will trigger a deceleration command 10ms in advance. This prediction mechanism makes up for the physical delay between sensor sampling and robotic arm execution, ensuring that momentum unloading is completed before the soft bag is irreversibly damaged.

[0078] Furthermore, in the impedance control stage, in order to find the optimal balance between "fast response" and "compliant protection", the system no longer uses a fixed mapping function, but adopts a real-time optimization algorithm based on gradient descent to dynamically calculate stiffness and damping coefficient. Construction of multi-objective cost function: The first term is the contact force stability index, which is calculated by the square of the difference between the real-time contact force and the expected contact force, and then multiplied by the force control weight coefficient (set to 0.7 to prioritize safety). The second item is the trajectory tracking accuracy index, which is calculated by multiplying the square of the difference between the actual position and the expected trajectory position by the position weight coefficient (set to 0.3 to take efficiency into account); the value of the cost function is the result of the weighted sum of the above two items. Furthermore, the calculation and updating of parameter gradients. Within each 5ms optimization cycle, the processor uses the numerical difference method to calculate the partial derivatives (i.e., gradients) of the cost function with respect to the stiffness coefficient and damping coefficient. The system sets the learning rate to 0.01; the update logic is: the new stiffness coefficient equals the current stiffness coefficient minus (learning rate multiplied by stiffness gradient value); the damping coefficient is updated similarly. To prevent parameter drift from causing system instability, the algorithm sets hard constraint boundaries: the stiffness coefficient is limited to between 500 N / m and 3000 N / m, and the damping coefficient is limited to between 50 Ns / m and 200 Ns / m. Furthermore, dynamic filtering and command generation specifically involve the optimal stiffness and damping parameters calculated by the optimization algorithm being immediately substituted into a dynamic filter constructed from a second-order impedance model. This filter processes the input desired velocity vector to generate an optimal position correction that balances minimizing contact force fluctuations and minimizing position errors. This method enables the stacker crane to finely adjust the stiffness of its arm in real time based on changes in the softness or hardness of the handle, just like a skilled worker when grabbing soft packages filled with liquid. This prevents the liquid from shaking violently and ensures the neatness of the stacking.

[0079] Although embodiments of the invention have been shown and described, it will be understood by those skilled in the art that various changes, modifications, substitutions and alterations can be made to these embodiments without departing from the principles and spirit of the invention, the scope of which is defined by the appended claims and their equivalents.

Claims

1. A multi-dimensional state observation based adaptive control system for industrial stacker cranes, characterized by, Comprise: State solving module: decoupling processing of the original observation data in the feature region, 6D pose solving for the rigid state subspace to obtain the position deviation matrix, and manifold geometry solving for the flexible state subspace to obtain the deformation state variable; Dynamic modeling module: monitoring the displacement response sequence at the top of the stacker column, online identifying the first-order natural frequency and damping ratio, and establishing a second-order impedance contact model containing virtual mass, damping coefficient and stiffness coefficient; Impedance control module: mapping the position deviation matrix into the expected speed vector of the servo system, adaptively adjusting the stiffness coefficient and damping coefficient of the second-order impedance contact model according to the deformation state variable, and generating the initial control instruction by using the adjusted second-order impedance contact model to modulate the expected speed vector; Estimation shaping module: constructing a Kalman state observer to perform state estimation on the initial control instruction and generate the modified speed control instruction; Based on the identified first-order natural frequency and damping ratio, a zero-vibration input shaper is constructed to perform time-domain convolution modulation on the modified speed control instruction to generate the shaping execution instruction.

2. The adaptive control system for industrial stacker based on multi-dimensional state observation according to claim 1, characterized in that, The acquisition of the original observation data comprises: Sending a time-sharing trigger signal to a dual-channel optical sensing unit, activating a first channel to project near-infrared speckle structured light, analyzing depth field information using the long-wave low-scattering characteristics of near-infrared light spectrum, then activating a second channel to project linearly polarized blue light, cooperating with an orthogonal polarization analyzer at the imaging end, collecting surface optical property information using the polarization extinction principle, converting the depth field information into a three-dimensional space coordinate model based on a pre-calibrated intrinsic parameter matrix and a hand-eye extrinsic parameter matrix, and correcting the distortion of the surface optical property information; performing spatial coordinate registration to map the surface optical property value to the corresponding three-dimensional space coordinate point to generate spatial state data as the original observation data.

3. The adaptive control system for industrial stacker based on multi-dimensional state observation according to claim 1, characterized in that, The state solving module is specifically configured to: Extract the three-dimensional space coordinate points and surface optical property values of the original observation data, establish a feature tensor, input it into a semantic segmentation neural network model, extract local topological features and surface attribute features, and calculate the state category probability of each data point; according to the maximum a posteriori probability principle, assign a state label to each data point to generate a rigid and flexible classification mask; apply the flexible classification mask to perform spatial filtering on the original observation data to separate the rigid state subspace and the flexible state subspace.

4. The adaptive control system for industrial stacker based on multi-dimensional state observation according to claim 1, characterized in that, The state solving module is specifically configured to: For the rigid state subspace, the iterative closest point registration algorithm is used to align the features of the spatial data in the region with the pre-set rigid standard model, calculate the rotation and translation transformation matrix of the current coordinate system relative to the ideal operation coordinate system, and obtain the position deviation matrix; for the flexible state subspace, the elliptic cylindrical surface manifold fitting algorithm is used to reconstruct the surface geometric topological structure of the region, calculate the long-short axis ratio and principal axis inclination parameters of the flexible operation object, and compare them with the pre-set object standard form to obtain the deformation state variable.

5. The adaptive control system for industrial stacker based on multi-dimensional state observation according to claim 1, characterized in that, The dynamic modeling module is specifically configured to: The top horizontal vibration displacement of the capturing column is captured to construct a time domain response sequence, a fast Fourier transform is performed to extract a first order natural frequency and a damping ratio is estimated by using a logarithmic attenuation method; a second order linear differential equation describing dynamic mapping of the position correction and the contact force is established as the impedance contact model, in which a virtual mass, a damping coefficient and a stiffness coefficient are respectively set to represent inertia response, energy dissipation rate and elastic restoring force of the system.

6. The adaptive control system for industrial stacker based on multi-dimensional state observation according to claim 1, characterized in that, The impedance control module is specifically configured to: A Jacobian matrix is constructed by using a position-based visual servoing control law, the position deviation matrix is taken as a feedback input, the expected speed vector of each servo axis is calculated through inverse kinematics solution and proportional gain adjustment, a nonlinear decay mapping function of the deformation state variable and the stiffness coefficient is established, the driving stiffness coefficient decreases along the mapping function curve when the deformation state variable increases, and the damping coefficient is adjusted synchronously according to the stiffness coefficient changing in real time based on a constant damping ratio constraint; The expected speed vector is taken as a reference input by using an admittance control structure, the adjusted stiffness coefficient and the damping coefficient are substituted into the second order impedance contact model to construct a dynamic filter, the expected speed vector is filtered to calculate the position and speed correction, and the initial control instruction is generated.

7. The adaptive control system for industrial stacker based on multi-dimensional state observation according to claim 1, characterized in that, The estimation shaping module is specifically configured to: A kinematics recursive model is constructed by using high-frequency sampling servo encoder feedback data, time update calculation of Kalman filtering is performed to obtain a prior state estimation value, the position deviation matrix with low-frequency sampling is taken as an observation vector to perform measurement update calculation of Kalman filtering, the prior state estimation value is corrected by using visual observation residual, and a corrected speed control instruction is generated, and a zero-vibration pulse sequence is constructed by using the identified first order natural frequency and damping ratio to calculate a pulse time interval and an amplitude coefficient; The zero-vibration pulse sequence is taken as a kernel function to perform time domain convolution operation on the corrected speed control instruction to obtain a shaping execution instruction containing a stepped pulse sequence.

8. The adaptive control system for industrial stacker based on multi-dimensional state observation according to claim 4, characterized in that, The impedance control module further includes virtual potential field guiding logic, which is specifically configured to: The surface geometry topology is received to construct a virtual repulsive potential field, and the maximum profile envelope surface of the flexible cargo is defined as the zero distance boundary of the potential field; a positive correlation mapping of the potential field intensity gradient and the deformation state variable is established, and the repulsive strength of the potential field is dynamically improved when the deformation state variable increases; in the process of generating the initial control instruction, the Euclidean distance gradient between the end effector and the potential field boundary is calculated in real time to generate a virtual repulsive speed vector and superimpose it on the expected speed vector, and a virtual reverse damping force is generated when the end effector approaches the cargo deformation boundary.

Images