Human–Robot Skill Transferring and Inverse Velocity Admittance Control for Soft Tissue Cutting Tasks

Abstract

Robotic meat cutting is increasingly in demand in meat industries due to safety issues, labor shortages, and inefficiencies. This paper proposes a multi-demonstration human–robot skill transfer framework to address the flexible and generalized cutting of sheep hindquarters with complex 3D anatomy structures by imitating humans. To improve the generalization with meat sizes and demonstrations and extract target cutting behaviors, multi-demonstrations of cutting are encoded into low-dimension latent space through principal components analysis (PCA), Gaussian mixture model (GMM), and Gaussian mixture regression (GMR). To improve the robotic cutting flexibility and the cutting behavior reproducing accuracy, this study combines a modified dynamic movement primitive (DMP) high-level behavior generator with the low-level joints admittance control (AC) through real-time inverse velocity (IV) kinematics solving and constructs the IVAC-DMP control module. The experimental results show that the maximum residual meat thickness in the sheep hindquarter cutting of sample 1 is 3.1 mm, and sample 2 is 3.8 mm. The residual rates of samples 1 and 2 are 5.6% and 4.8%. Both meet the requirements for sheep hindquarter separation. The proposed framework is advantageous for harvesting high-value meat products and providing a reference technique for robot skill learning in interaction tasks.

1. Introduction

Meat cutting is a crucial component of the meat industry’s value chain. This process involves breaking down a carcass into primal parts and precisely harvesting meat cuts with varying values [1,2]. Skilled workers with a profound understanding of animal anatomy are vital to guaranteeing the final cuts’ consistency and quality and determining their market prices. However, the traditional manual cutting methods have several issues, including biological variation, training, sanitation concerns, and harsh working conditions [3]. With the recent labor shortages, supply inefficiencies caused by COVID-19, and other limitations, the meat industry has seen a considerable increase in the need for automation. Robotic solutions offer a promising approach to address these challenges by enabling human-like cutting operations, reducing labor costs, improving food safety, and promoting a sustainable meat supply [4].

Robotic cutting is achieved through three main methods [5], i.e., primal sawing, 2D-space, and 3D-space cutting. The primal sawing method is efficient in primal cutting and middle pieces separation [6]. However, the simple sawing actions are insufficient for harvesting meat cuts from forequarters, hindquarters, rib cages, and carcass flanks. The 2D-space methods have been primarily adopted for precisely cutting primal pieces [7,8], such as carcass flanks or rib cages, utilizing deep learning and multi-DoF robot systems. The 2D-space methods are highly limited in cutting meat with more complex 3D anatomy structures, e.g., sheep hindquarters. The motion planning and control of robot cutting for sheep hindquarter separation is difficult due to several factors, i.e., cutting scenario variability, cutting close to the bone, variable bone size, and complex anatomic structures [9]. Compliant cutting in a 3D space can precisely separate meat tissue from bones or different meat tissues from each other. Xie et al. [10] proposed an admittance control method to cut along the ischial bone with comparable efficiency and accuracy to humans for sheep hindquarter cutting. However, the cutting path was constrained on several planes around the bone, which limited the flexibility of the cutting. Meat and Livestock Australia (MLA) proposed a force-controlled robot for sheep hindquarter cutting [11]. The robot equipped a blade with an optimized shape to lower the cutting forces. However, the generalization of the cutting motions is not addressed. Furthermore, a robotic prototype was proposed to cut the shoulders of the sheep forequarter [12] through an oscillating cutting motion and force control. It is a problem to generalize the oscillating cutting pattern between the forequarters and hindquarters. In addition, Nabil et al. [13] utilized a finite element model to accurately predict the deformation of the meat tissue during cutting and control the robot’s motion. Although the study allows for precise tracking of the fascial tissues between the meat tissues of a beef hind leg, the bones were not considered. More work is needed to achieve flexible and generalized motion planning while coupling with the force control.

Robots have been trained to execute flexible motions in interaction tasks based on the human–robot skill transfer (HRST) framework [14]. The framework involves a high-level nonlinear dynamical system such as a dynamic movement primitive (DMP) and a low-level compliant controller [15]. The nonlinear dynamical system combines a set of radial basis functions (RBFs) that mimic the muscle unit behaviors in humans to generate flexible motion trajectories [16]. Meanwhile, the low-level controller is an impedance controller(IC) or admittance controller (AC) for the robot joints, which handles the interaction between the robot and objects [17,18]. The DMP model has been coupled with feedback terms related to the interaction force to form a Cartesian space AC [19,20] in grinding, handling, and robot collaboration tasks. However, the method was highly limited by inverse kinematics solving (IKS), since the desired trajectories were modulated online. In the primary research [21], we utilized the DMP model to generate the desired movement in sheep hindquarter separation and the inverse velocity (IV) algorithm for online IKS of the robot with AC. However, the low-level joint control has no feedback to the DMP model, which limits the flexibility of the method. The studies mentioned above consider joint space control and high-level planning separately, so the robotic cutting flexibility and the cutting behavior reproducing accuracy are highly limited. More studies combining low-level joint AC and the high-level DMP behavior system by real-time IV kinematics solving are in demand to improve the integration, consistency, and efficiency of the collaboration between the two systems.

Moreover, The studies mentioned above are single-demonstrated, whose generalization is a crucial issue due to the cutting objects’ variable types, locations, and dimensions. The generalization can be addressed through multi-demonstration HRST. Gams et al. [22] proposed the compliant movement primitive (CMP) model to learn the interactive forces and motion trajectories and utilized the Gaussian process regression (GPR) model to encode multi-demonstrations. Low motion tracking error and compliance control were achieved. In addition, Wu and Billard [23] combined the Stable Estimator of Dynamical Systems) and Gaussian mixture model (GMM) and Gaussian mixture regression (GMR) to learn a compliant cutting behavior model from multiple demonstrations. The model was then integrated into the IC of the robot for cutting single-layered skin-like soft tissue and responding to human interactions. The studies mentioned above are implemented under an IC framework. An expensive lightweight backdrivable robot arm is needed to produce joint torques precisely. In the task of cutting hindquarters, HRST through multi-demonstration on a robot with AC and real-time IV remains an open problem.

To address the efficient collaboration between human-like behavior generation and joint-space compliant control and achieve flexible and generalized hindquarter cutting, this paper proposes a multi-demonstration HRST framework. To improve the generalization with meat sizes and demonstrations and extract target cutting behaviors, multi demonstrations of cutting are encoded into a low-dimension latent space through principal components analysis (PCA), Gaussian mixture model (GMM), and Gaussian mixture regression (GMR). To improve the robotic cutting flexibility and the cutting behavior reproducing accuracy, this study combines a modified dynamic movement primitive (DMP) high-level behavior generator with the low-level joints admittance control (AC) through real-time inverse velocity (IV) kinematics solving and constructs the IVAC-DMP control module. The proposed framework is evaluated in four cutting scenarios. The proposed method has the potential to reduce the cutting force, minimize meat loss, and increase meat yield, which is advantageous for harvesting high-value meat products. Additionally, it can serve as a reference technique for flexible and generalized robot skill learning in interaction tasks.

The rest of the paper is structured as follows. The proposed framework is presented in Section 2. Section 2.2 presents the data acquisition and preprocessing. The multi-demonstration human–robot skill transfer with GMM-GMR is then explained in Section 2.3. Moreover, Section 2.4 demonstrates the IVAC of robot joints. The evaluation and the discussion are presented in Section 3. The conclusion is provided in Section 4.

2. Materials and Methods

This section provides a detailed description of the proposed framework for imitating and controlling the cutting of soft tissue by robots.

2.1. System Overview

2.1.1. Cutting Behavior Representation

The cutting behavior is represented by ${\mathbf{h}}_{\mathrm{d}}\left(t\right)$,

$$\begin{array}{cc}\hfill {\mathbf{H}}_{\mathrm{d}}& =[{\mathbf{h}}_{\mathrm{d}}\left(t_1\right),{\mathbf{h}}_{\mathrm{d}}\left(t_2\right),\cdots ,{\mathbf{h}}_{\mathrm{d}}\left(t_M\right)],\hfill \\ \hfill {\mathbf{F}}_{\mathrm{d}}& =[{\mathbf{f}}_{\mathrm{d}}\left(t_1\right),{\mathbf{f}}_{\mathrm{d}}\left(t_2\right),\cdots ,{\mathbf{f}}_{\mathrm{d}}\left(t_M\right)],\hfill \\ \hfill {\mathbf{T}}_{\mathrm{d}}& =[t_1,t_2,\cdots ,t_M],\hfill \\ \hfill {\mathbf{h}}_{\mathrm{d}}\left(t\right)& ={[{\mathbf{x}}_{\mathrm{d}}\left(t\right),{\mathbf{q}}_{\mathrm{d}}\left(t\right)]}^{\mathrm{T}},\hfill \\ \hfill {\mathbf{x}}_{\mathrm{d}}\left(t\right)& ={[x_{\mathrm{d},\mathrm{x}}\left(t\right),x_{\mathrm{d},\mathrm{y}}\left(t\right),x_{\mathrm{d},\mathrm{z}}\left(t\right)]}^{\mathrm{T}},\hfill \\ \hfill {\mathbf{q}}_{\mathrm{d}}\left(t\right)& ={[q_{\mathrm{d},\mathrm{w}}\left(t\right),q_{\mathrm{d},\mathrm{x}}\left(t\right),q_{\mathrm{d},\mathrm{y}}\left(t\right),q_{\mathrm{d},\mathrm{z}}\left(t\right)]}^{\mathrm{T}},\hfill \\ \hfill {\mathbf{f}}_{\mathrm{d}}\left(t\right)& ={[f_{\mathrm{d},\mathrm{x}}\left(t\right),f_{\mathrm{d},\mathrm{y}}\left(t\right),f_{\mathrm{d},\mathrm{z}}\left(t\right),{\tau }_{\mathrm{d},\mathrm{x}}\left(t\right),{\tau }_{\mathrm{d},\mathrm{x}}\left(t\right),{\tau }_{\mathrm{d},\mathrm{x}}\left(t\right)]}^{\mathrm{T}},\hfill \end{array}$$

(1)

where t is the time, and ${\mathbf{H}}_{\mathrm{d}}$ is the collected time series of cutting behavior. M is the number of time points in the time series. ${\mathbf{F}}_{\mathrm{d}}$ is the collected time series of cutting force. ${\mathbf{T}}_{\mathrm{d}}$ is the time points series. The position of the cutting tool in the world coordinate system is represented by ${\mathbf{x}}_{\mathrm{d}}\left(t\right)$, which has $x_{\mathrm{d},\mathrm{x}}\left(t\right),x_{\mathrm{d},\mathrm{y}}\left(t\right),x_{\mathrm{d},\mathrm{z}}\left(t\right)$ as coordinates for the $\mathrm{x},\mathrm{y},\mathrm{z}$ axes, respectively. The quaternion of the cutting tool orientation in the world coordinate system is represented by ${\mathbf{q}}_{\mathrm{d}}\left(t\right)$, which consists of $q_{\mathrm{d},\mathrm{w}}\left(t\right)$ as the scaler part and $q_{\mathrm{d},\mathrm{x}}\left(t\right),q_{\mathrm{d},\mathrm{y}}\left(t\right),q_{\mathrm{d},\mathrm{z}}\left(t\right)$ as the imaginary part. The cutting force/torque measured by the F/T sensor is represented by ${\mathbf{f}}_{\mathrm{d}}\left(t\right)$, where $f_{\mathrm{d},\mathrm{x}}\left(t\right),f_{\mathrm{d},\mathrm{y}}\left(t\right),f_{\mathrm{d},\mathrm{z}}\left(t\right)$ are the decoupled measured forces of $\mathrm{x},\mathrm{y},\mathrm{z}$ while ${\tau }_{\mathrm{d},\mathrm{x}}\left(t\right),{\tau }_{\mathrm{d},\mathrm{x}}\left(t\right),{\tau }_{\mathrm{d},\mathrm{x}}\left(t\right)$ are the decoupled measured torques of the $\mathrm{x},\mathrm{y},\mathrm{z}$ axes, respectively.

2.1.2. Process of Training and Control

The process of imitation learning and control is divided into three stages, as shown in Figure 1. By giving multiple groups of human behavior time-series samples in the same cutting scenario, we train a GMM in the probability distribution of human behavior at each time point. Then, we generate the target cutting behavior through GMR with the trained GMM. The target cutting behavior is transformed into a DMP system in this work to generate the target cutting behavior (position, orientation, and feedforward cutting force) of the robot end effector in the base coordinate system. Finally, we input the target cutting behavior to the IVAC module to control the joint angles of the robot during cutting to achieve human-like flexible cutting.

2.1.3. Experiment Platform

As is shown in Figure 2, the experiment platform includes a loading device, a cutting tool, a visual motion capture system, and a multi-degree-of-freedom robot arm. The cutting tool is constructed by connecting a meat separation blade with a six-axis F/T sensor using a metal frame. The ArUco augmented reality marker [24] is fixed to the left side of the cutting tool, and two intel RealSense D435i cameras track its pose to form a visual motion capture system. The D435i cameras detect the ArUco marker at a frequency of 30 Hz. They verify the marker’s encoding ID and utilize the corner detection tool in OpenCV to determine the marker’s position and orientation in the camera’s coordinate system. This process enables the camera to track the marker accurately. The marker’s ID is 123. This unique ID distinguishes and verifies other QR codes within view. The marker is 60 mm × 60 mm, ensuring clear visibility from a preset distance. The marker is a finely machined anodized aluminum sheet with distinct corners and low reflection. Meanwhile, the F/T sensor measures the decoupled forces and torques in the x, y, and z axes.

The F/T sensor has flanges at the front and rear; the front flange is connected to the metal frame of the cutting tool, and the rear flange can be connected with the flange of the robot end or with a 3D-printed grip. When the rear flange connects with the grip, the demonstrator holds the cutting tool by the grip to cut. Otherwise, the robot arm manipulates the cutting tool to cut.

2.2. Data Acquisition and Preprocessing

2.2.1. Data Acquisition

This study collects the cutting behavior and force data in four demonstrating scenarios to construct the datasets. The experiment was conducted in the College of Engineering, China Agricultural University with indoor white light and at a room temperature of 17 °C. The experiment materials include the sheep meat from the local (Beijing) meat plants and foam blocks. The raw data are calibrated, aligned, and then utilized to train the GMM-GMR model. To ensure the robot learns the exact human behaviors, this study only uses the actual demonstrations from real humans, and no data augmentation methods are employed. The datasets include 9570 data points (2348 from scenario 1, 884 from scenario 2, 1358 from scenario 3, and 4980 from scenario 4) from 38 demonstration samples (4 from scenario 1, 2 from scenario 2, 2 from scenario 3, and 30 from scenario 4). Scenario 1 is cutting the medial hind-leg muscles attached with the ischium ventral; 2 is cutting the lateral hind-leg muscles attached with the ischium dorsal; 3 is cutting the lateral hind-leg muscles attached with the ilium dorsal; 4 is to segment a foam block pasted on a rigid surface.

2.2.2. Data Calibration and Filtering

The raw data are calibrated to obtain the transform matrices between the robot base coordinate system W, the robot end mounting flange center coordinate system E, the F/T sensor coordinate system S, the marker coordinate system T, the Camera 1 coordinate system $C_1$, and the Camera 2 coordinate system $C_2$. ${\mathbf{T}}_{WE}$ represents the transform matrix of W and E. ${\mathbf{T}}_{ES}$ represents the transform matrix of E and S. ${\mathbf{T}}_{C_{1}T}$ represents the transform matrix of $C_1$ and T. ${\mathbf{T}}_{C_{2}T}$ represents the transform matrix of $C_2$ and T. ${\mathbf{T}}_{W{C}_1}$ represents the transform matrix of W and $C_1$. ${\mathbf{T}}_{W{C}_2}$ represents the transform matrix of W and $C_2$. ${\mathbf{T}}_{TE}$ represents the transform matrix of T and E. The data calibration process includes seven steps:

• ${\mathbf{T}}_{WE}$ is the homogeneous matrix representation of E’s pose in W. The pose is read from the law-level embedded controller of the robot. Since the end mounting flange and the F/T sensor are connected by a mechanical structure, ${\mathbf{T}}_{ES}$ can be pre-obtained.
• ${\mathbf{T}}_{C_{1}T}$ is the homogeneous matrix representation of T’s pose in $C_1$, and ${\mathbf{T}}_{C_{2}T}$ is the homogeneous matrix representation of T’s pose in $C_2$. ${\mathbf{T}}_{C_{1}T}$ and ${\mathbf{T}}_{C_{2}T}$ are provided by the ArUco marker detection algorithm. Furthermore, ${\mathbf{T}}_{WE}$, ${\mathbf{T}}_{C_{1}T}$, and ${\mathbf{T}}_{C_{2}T}$, ${\mathbf{T}}_{W{C}_1}$, ${\mathbf{T}}_{W{C}_2}$, ${\mathbf{T}}_{T{E}_1}$, and ${\mathbf{T}}_{T{E}_2}$ are obtained through hand–eye calibration.
• We convert the measured pose of T on $C_1$ and $C_2$ to the pose of E on W with Equation (2) using the calibrated ${\mathbf{T}}_{W{C}_1}$, ${\mathbf{T}}_{W{C}_2}$, ${\mathbf{T}}_{T{E}_1}$, and ${\mathbf{T}}_{T{E}_2}$,

  $$\begin{array}{cc}\hfill {\mathbf{T}}_{W{E}_1}& ={\mathbf{T}}_{W{C}_1}{\mathbf{T}}_{C_{1}T}{\mathbf{T}}_{T{E}_1},\hfill \\ \hfill {\mathbf{T}}_{W{E}_2}& ={\mathbf{T}}_{W{C}_2}{\mathbf{T}}_{C_{2}T}{\mathbf{T}}_{T{E}_2}.\hfill \end{array}$$

  (2)
• We rewrite the homogeneous matrices ${\mathbf{T}}_{W{E}_1}$ and ${\mathbf{T}}_{W{E}_2}$ into pose vectors and obtain the position ${\mathbf{x}}_{\mathrm{d}}$ and quaternion ${\mathbf{q}}_{\mathrm{d}}$ of the cutting tool in the world coordinate system by Kalman filtering the pose vectors from the two camera sources.
• We convert the cutting forces measured on the coordinate system S to the coordinate system E by (3) using the calibrated ${\mathbf{T}}_{SE}$,

  $${\mathbf{f}}_{\mathrm{e}}={\left[{\mathbf{Ad}}_{{\mathbf{T}}_{SE}}\right]}^{\mathrm{T}}{\mathbf{f}}_{\mathrm{s}},$$

  (3)

  where ${\mathbf{f}}_{\mathrm{e}}$ denotes the cutting force on E, ${\mathbf{f}}_{\mathrm{s}}$ denotes the cutting force on S, and $\left[{\mathbf{Ad}}_{{\mathbf{T}}_{SE}}\right]$ denotes the adjoint representation of ${\mathbf{T}}_{SE}$.
• To ensure accurate cutting force measurements, we use a gravity-compensated calibration method [25] to eliminate the impact of the cutting tool’s self-weight on the collected data. The resulting cutting force is denoted as ${\mathbf{f}}_{\mathrm{d}}$. The components of this force, $f_{\mathrm{d},\mathrm{x}}$, $f_{\mathrm{d},\mathrm{y}}$, and $f_{\mathrm{d},\mathrm{z}}$, represent different directions. $f_{\mathrm{d},\mathrm{x}}$ points towards the side of the blade, $f_{\mathrm{d},\mathrm{z}}$ points towards the tip, and $f_{\mathrm{d},\mathrm{y}}$ points towards the edge, all within the blade plane.
• Dynamic time warping (DTW) is utilized to align the sample time-series length in the same scenario in a Python environment.

2.2.3. Latent Space Representation of Dataset

In order to improve learning efficiency, this study utilizes the encoding of the high-dimension time series in the dataset into the one-dimension latent variable series ${\mathbf{Z}}_{\mathrm{cat}}$. First, the samples in the dataset are concatenated, as shown in Equation (4),

$$\begin{array}{cc}\hfill {\mathbf{H}}_{\mathrm{cat}}& =[{\mathbf{H}}_{\mathrm{d},1},{\mathbf{H}}_{\mathrm{d},2},\dots ,{\mathbf{H}}_{\mathrm{d},L}],\hfill \\ \hfill {\mathbf{T}}_{\mathrm{cat}}& =[{\mathbf{T}}_{\mathrm{d},1},{\mathbf{T}}_{\mathrm{d},2},\dots ,{\mathbf{T}}_{\mathrm{d},L}],\hfill \end{array}$$

(4)

where ${\mathbf{H}}_{\mathrm{cat}}$ denotes the concatenated series of cutting behavior. ${\mathbf{T}}_{\mathrm{cat}}$ denotes the concatenated series of time points. L denotes the number of samples. ${\mathbf{H}}_{\mathrm{d},1},{\mathbf{H}}_{\mathrm{d},2},\dots ,{\mathbf{H}}_{\mathrm{d},L}$ are the samples. Second, we perform a principal components analysis (PCA) on ${\mathbf{H}}_{\mathrm{cat}}$. The PCA extracts the correlation and importance information among different axes of the cutting tool in the latent variables. The feature vector with the largest eigenvalue obtained by PCA is as follows,

$${\mathbf{W}}_{\mathrm{pca}}=\left[w_{\mathrm{pca},1},w_{\mathrm{pca},2},\dots ,w_{\mathrm{pca},6}\right],$$

(5)

where ${\mathbf{W}}_{\mathrm{pca}}$ denotes the feature vector with elements $w_{\mathrm{pca},1},w_{\mathrm{pca},2},\dots ,w_{\mathrm{pca},6}$. Then, ${\mathbf{H}}_{\mathrm{cat}}$ is transformed into its latent variable representation using Equation (6),

$${\mathbf{Z}}_{\mathrm{cat}}={\mathbf{W}}_{\mathrm{pca}}{\mathbf{H}}_{\mathrm{cat}},$$

(6)

where ${\mathbf{Z}}_{\mathrm{cat}}$ denotes the latent variable series with a length of $L\times M$.

2.3. Learning of Cutting Behaviors with GMM and GMR

2.3.1. Learning from Multi-demonstrations of Cutting Behaviors with GMM

The GMM is a statistical model for representing a given population as several Gaussian-distributed subpopulations. Therefore, it is suitable for learning cutting behavior with multi-demonstrations from the same scenario [26]. In this study, the GMM is trained on the dataset $\left\{{\mathbf{Z}}_{\mathrm{cat}},{\mathbf{T}}_{\mathrm{cat}}\right\}$, with the state vector consisting of the cutting behavior latent variable $z_{\mathrm{cat}}$ and time $t_{\mathrm{cat}}$, as shown in Equation (7),

$$\begin{array}{cc}\hfill p\left(\mathit{\xi }\right)& \sim \sum _{k=1}^K{\pi }_k{\mathcal{N}}_k\left({\mathit{\mu }}_k,{\mathbf{\Sigma }}_k\right),\mathrm{with}\sum _{k=1}^K{\pi }_k=1,\hfill \\ \hfill \mathit{\xi }& =\left[\begin{array}{c}{t}_{\mathrm{cat}}\\ z_{\mathrm{cat}}\end{array}\right],{\mathit{\mu }}_k=\left[\begin{array}{c}{\mu }_{\mathrm{t},k}\\ {\mu }_{\mathrm{z},k}\end{array}\right],{\mathbf{\Sigma }}_k=\left[\begin{array}{cc}{\mathrm{\Sigma }}_{\mathrm{t},k}& {\mathrm{\Sigma }}_{\mathrm{tz},k}\\ {\mathrm{\Sigma }}_{\mathrm{zt},k}& {\mathrm{\Sigma }}_{\mathrm{z},k}\end{array}\right],\hfill \end{array}$$

(7)

where $p\left(\mathit{\xi }\right)$ is the probability distribution of state vectors, and $p\left(\mathit{\xi }\right)$ is the mixture of K multivariant Gaussians. ${\mathcal{N}}_k\left({\mathit{\mu }}_k,{\mathbf{\Sigma }}_k\right)$ represents the kth Gaussian with a center ${\mathit{\mu }}_k$, a covariance matrix ${\mathbf{\Sigma }}_k$, and a weight ${\pi }_k$.

The values of ${\pi }_k$, ${\mathit{\mu }}_k$, and ${\mathbf{\Sigma }}_k$ are obtained by training on input data through the expectation maximization (EM) method. The EM algorithm iteratively improves the similarity between the GMM and the realistic distribution of data samples to achieve maximum likelihood estimation of the sample distribution. The iterative optimization process involves initialization, E-step, and M-step:

• The mixture model consists of 10 Gaussians. To reduce the sensitivity of the EM algorithm to the selection of the initial values, we use k-means to initialize the centers and covariance matrices of the Gaussians. After initializing, we train the GMM with the dataset of the four cutting scenarios that were collected;
• The E-step calculates the intermediate variable ${\gamma }_{m,k}$ using current ${\mu }_k$ and ${\mathbf{\Sigma }}_k$, as is shown in Equation (8),

  $${\gamma }_{m,k}=\frac{{\pi }_k{\mathcal{N}}_k\left({\mathit{\xi }}_m|{\mathit{\mu }}_k,{\mathbf{\Sigma }}_k\right)}{{\sum }_{k=1}^K{\pi }_k{\mathcal{N}}_k\left({\mathit{\xi }}_m|{\mathit{\mu }}_k,{\mathbf{\Sigma }}_k\right)}.$$

  (8)

  where ${\gamma }_{m,k}$ denotes the condition probability of ${\mathit{\xi }}_m$ generating from the kth Gaussian;
• Then, the M-step updates ${\pi }_k$, ${\mathit{\mu }}_k$, and ${\mathbf{\Sigma }}_k$ using ${\gamma }_{m,k}$ calculated in the E-step, as is shown in Equation (9),

  $$\begin{array}{cc}\hfill {\pi }_k^{\mathrm{new}}& =\frac{{\sum }_{m=1}^M{\gamma }_{m,k}}{M},\hfill \\ \hfill {\mathit{\mu }}_k^{\mathrm{new}}& =\frac{{\sum }_{m=1}^M{\gamma }_{m,k}{\mathit{\xi }}_m}{{\sum }_{m=1}^M{\gamma }_{m,k}},\hfill \\ \hfill {\mathbf{\Sigma }}_k^{\mathrm{new}}& =\frac{{\sum }_{m=1}^M{\gamma }_{m,k}\left({\mathit{\xi }}_m-{\mathit{\mu }}_k\right){\left({\mathit{\xi }}_m-{\mathit{\mu }}_k\right)}^{\mathrm{T}}}{{\sum }_{m=1}^M{\gamma }_{m,k}},\hfill \end{array}$$

  (9)

  where ${\pi }_k^{\mathrm{new}}$ denotes the weight after updating. ${\mathit{\mu }}_k^{\mathrm{new}}$ denotes the center after updating. ${\mathbf{\Sigma }}_k^{\mathrm{new}}$ denotes the covariance matrix after updating.

2.3.2. Robotic Cutting Behavior Generation with GMR

As is shown in Figure 3, this study employs GMR to generate the target latent space trajectory ${\mathbf{Z}}_{\mathrm{ref}}$ for the cutting behaviors. In Figure 3a–d, the green-colored dots represent the points in the datasets, and thick red lines represent regression trajectories generated by the GMR. The ellipses represent the Gaussians in the GMM. The radius and angle of every ellipse represent a Gaussian covariance matrix ${\mathbf{\Sigma }}_k$, and the crosses represent the mean values of the Gaussians. Figure 3e–h show the corresponding cutting paths. The latent variable trajectories are almost averaged out by the GMM-GMR model in the four scenarios. Significant differences are found in the GMM and GMR trajectories learned in the four scenarios, since the cutting paths vary.

The regression function is calculated using the following equation:

$$z_{\mathrm{ref}}\left(t\right)={\widehat{\mu }}_{\mathrm{z}}\left(t\right)=\sum _{k=1}^K{w}_k\left(t\right){\widehat{\mu }}_{\mathrm{z},k}\left(t\right),$$

(10)

with,

$$\begin{array}{cc}\hfill {\widehat{\mu }}_{\mathrm{z},k}\left(t\right)& ={\mu }_{\mathrm{z},k}+{\Sigma }_{\mathrm{z},k}{\Sigma }_{\mathrm{t},k}^{-1}(t-{\mu }_{\mathrm{t},k}),\hfill \\ \hfill w_k\left(t\right)& =\frac{{\pi }_k\mathcal{N}\left(t\right|{\mu }_{\mathrm{t},k},{\Sigma }_{\mathrm{t},k})}{{\sum }_{k=1}^K{\pi }_k\mathcal{N}\left(t\right|{\mu }_{\mathrm{t},k},{\Sigma }_{\mathrm{t},k})},\hfill \end{array}$$

(11)

where $z_{\mathrm{ref}}\left(t\right)$ denotes the regression value of the latent variable at time t. ${\widehat{\mu }}_{\mathrm{z},k}\left(t\right)$ denotes the regression value of ${\mu }_{\mathrm{z},k}$ at time t, and $w_k$ denotes the weight. ${\mu }_{\mathrm{z},k}$, ${\mathrm{\Sigma }}_{\mathrm{t},k}$, and ${\mathrm{\Sigma }}_{\mathrm{z},k}$ are parameters defined in Equation (7). In this study, the values of the latent variables were calculated using Formula (10) at each time point, and the reference cutting behavior time series ${\mathbf{Z}}_{\mathrm{ref}}$ in latent space was obtained,

$${\mathbf{Z}}_{\mathrm{ref}}=[z_{\mathrm{ref}}\left(t_1\right),z_{\mathrm{ref}}\left(t_2\right),\cdots ,z_{\mathrm{ref}}\left(t_M\right)],$$

(12)

2.3.3. Post-Processing

The latent variable series of cutting behavior obtained by GMR were post-processed to generate the reference cutting behavior time series. First, by decoding the sequence of latent variables in (12), the reference cutting behavior (positions, quaternions, and forces) time series is generated,

$$\begin{array}{cc}\hfill {\mathbf{H}}_{\mathrm{ref}}& =[h_{\mathrm{ref}}\left(t_1\right),h_{\mathrm{ref}}\left(t_2\right),\cdots ,h_{\mathrm{ref}}\left(t_M\right)],\hfill \\ \hfill {\mathbf{H}}_{\mathrm{ref}}& ={\mathbf{W}}_{\mathrm{pca}}^{\mathrm{T}}·{\mathbf{Z}}_{\mathrm{ref}},\hfill \\ \hfill h_{\mathrm{ref}}\left(t_i\right)& ={\mathbf{W}}_{\mathrm{pca}}^{\mathrm{T}}·z_{\mathrm{ref}}\left(t_i\right),\mathrm{with}i=1,\cdots ,M,\hfill \end{array}$$

(13)

where the symbol ${\mathbf{H}}_{\mathrm{ref}}$ refers to the time series of reference cutting behaviors. Each element $h_{\mathrm{ref}}\left(t_i\right)$ is a cutting behavior vector in the form of Equation (1) at time $t_i$. The symbol ${\mathbf{W}}_{\mathrm{pca}}$ refers to the eigenvector in Equation (5). The symbol $z_{\mathrm{ref}}\left(t_i\right)$ is a latent variable at time $t_i$. The time series has a length of M.

According to Equation (1), the time series ${\mathbf{H}}_{\mathrm{ref}}$ of reference cutting behaviors is divided into ${\mathbf{X}}_{\mathrm{ref}}$, ${\mathbf{Q}}_{\mathrm{ref}}$, and ${\mathbf{F}}_{\mathrm{ref}}$,

$$\begin{array}{cc}\hfill {\mathbf{X}}_{\mathrm{ref}}& =[{\mathbf{x}}_{\mathrm{ref}}\left(t_1\right),{\mathbf{x}}_{\mathrm{ref}}\left(t_2\right),\cdots ,{\mathbf{x}}_{\mathrm{ref}}\left(t_M\right)],\hfill \\ \hfill {\mathbf{Q}}_{\mathrm{ref}}& =[{\mathbf{q}}_{\mathrm{ref}}\left(t_1\right),{\mathbf{q}}_{\mathrm{ref}}\left(t_2\right),\cdots ,{\mathbf{q}}_{\mathrm{ref}}\left(t_M\right)],\hfill \\ \hfill {\mathbf{F}}_{\mathrm{ref}}& =[{\mathbf{f}}_{\mathrm{ref}}\left(t_1\right),{\mathbf{f}}_{\mathrm{ref}}\left(t_2\right),\cdots ,{\mathbf{f}}_{\mathrm{ref}}\left(t_M\right)],\hfill \end{array}$$

(14)

where ${\mathbf{X}}_{\mathrm{ref}}$ denotes the reference tool position sequence, ${\mathbf{Q}}_{\mathrm{ref}}$ denotes the reference tool orientation sequence, and ${\mathbf{F}}_{\mathrm{ref}}$ denotes the reference cutting force sequence.

Furthermore, the first- and second-order derivatives of the reference positions and quaternions are calculated,

$$\begin{array}{cc}\hfill {\mathbf{V}}_{\mathrm{ref}}& =[{\mathbf{v}}_{\mathrm{ref}}\left(t_1\right),{\mathbf{v}}_{\mathrm{ref}}\left(t_2\right),\cdots ,{\mathbf{v}}_{\mathrm{ref}}\left(t_M\right)],\hfill \\ \hfill {\dot{\mathbf{V}}}_{\mathrm{ref}}& =[{\dot{\mathbf{v}}}_{\mathrm{ref}}\left(t_1\right),{\dot{\mathbf{v}}}_{\mathrm{ref}}\left(t_2\right),\cdots ,{\dot{\mathbf{v}}}_{\mathrm{ref}}\left(t_M\right)],\hfill \\ \hfill {\Omega }_{\mathrm{ref}}& =[{\mathit{\omega }}_{\mathrm{ref}}\left(t_1\right),{\mathit{\omega }}_{\mathrm{ref}}\left(t_2\right),\cdots ,{\mathit{\omega }}_{\mathrm{ref}}\left(t_M\right)],\hfill \\ \hfill {\dot{\Omega }}_{\mathrm{ref}}& =[{\dot{\mathit{\omega }}}_{\mathrm{ref}}\left(t_1\right),{\dot{\mathit{\omega }}}_{\mathrm{ref}}\left(t_2\right),\cdots ,{\dot{\mathit{\omega }}}_{\mathrm{ref}}\left(t_M\right)],\hfill \end{array}$$

(15)

with

$$\begin{array}{cc}\hfill {\mathbf{v}}_{\mathrm{ref}}\left(t_i\right)& =\frac{{\mathbf{x}}_{\mathrm{ref}}\left(t_i\right)-{\mathbf{x}}_{\mathrm{ref}}\left(t_{i-1}\right)}{t_i-t_{i-1}},i>1,\hfill \\ \hfill {\dot{\mathbf{v}}}_{\mathrm{ref}}\left(t_i\right)& =\frac{{\mathbf{v}}_{\mathrm{ref}}\left(t_i\right)-{\mathbf{v}}_{\mathrm{ref}}\left(t_{i-1}\right)}{t_i-t_{i-1}},i>1,\hfill \\ \hfill {\mathit{\omega }}_{\mathrm{ref}}\left(t_i\right)& =\frac{{\mathbf{d}}_{\mathrm{q}}({\mathbf{q}}_{\mathrm{ref}}\left(t_{i-1}\right),{\mathbf{q}}_{\mathrm{ref}}\left(t_i\right))}{t_i-t_{i-1}},i>1,\hfill \\ \hfill {\dot{\mathit{\omega }}}_{\mathrm{ref}}\left(t_i\right)& =\frac{{\mathit{\omega }}_{\mathrm{ref}}\left(t_i\right)-{\mathit{\omega }}_{\mathrm{ref}}\left(t_{i-1}\right)}{t_i-t_{i-1}},i>1,\hfill \end{array}$$

(16)

where $t_i$ denotes the time with $t_1=0$, ${\mathbf{v}}_{\mathrm{ref}}\left(t_1\right)=0\mathrm{m}/\mathrm{s},{\dot{\mathbf{v}}}_{\mathrm{ref}}\left(t_1\right)=0\mathrm{m}/{\mathrm{s}}^2,{\mathit{\omega }}_{\mathrm{ref}}\left(t_1\right)=0\mathrm{rad}/\mathrm{s},\mathrm{and}{\dot{\mathit{\omega }}}_{\mathrm{ref}}\left(t_1\right)=0\mathrm{rad}/{\mathrm{s}}^2$. The unit quaternion representation is adopted to avoid the orientation singularity. ${\mathbf{d}}_{\mathrm{q}}({\mathbf{q}}_{\mathrm{ref}}\left(t_{i-1}\right),{\mathbf{q}}_{\mathrm{ref}}\left(t_i\right))$ denotes the orientation offset, which maps the quaternion offset into angular velocity space for derivation.

2.4. Admittance Control in Joint Space Coupling with DMPs

2.4.1. Dynamic Motion Primitives of Target Cutting Behaviors

(1) Activation Function: Like original DMPs, we convert time variables into phase space activation signals on $\left[0,1\right)$ through an activation function. We take the Sigmoid system (17) with upper and lower limits as the activation function

$$\tau \dot{s}={\alpha }_{s}s(1-s),$$

(17)

where s is the phase space activation signal, and ${\alpha }_s=2log(0.01/0.99)$. s starts at 0.99 and decreases to 0.01 at $t=\tau$ to simulate neural signal activation and attenuation.

(2) Modified DMPs model: first, considering the transformation relationship between quaternions and angular velocities, we extend the model proposed by Hoffmann et al. [27] to the quaternion space and build the following improved DMP model (18) for reference cutting behavior,

$${\tau }^2\left[\begin{array}{c}{\mathbf{a}}_{\mathrm{d}}\\ {\mathit{\alpha }}_{\mathrm{d}}\end{array}\right]=\left[\begin{array}{cc}{\mathbf{K}}_{\mathrm{p}}& \\ & {\mathbf{K}}_{\mathrm{q}}\end{array}\right]\left[\begin{array}{c}s\left({\mathbf{w}}_{\mathrm{p}}^{\mathrm{T}}\mathbf{g}-\mathbf{x}+{\mathbf{x}}_0\right)+\left(1-s\right)({\mathbf{x}}_g-\mathbf{x})\\ s\left({\mathbf{w}}_{\mathrm{q}}^{\mathrm{T}}\mathbf{g}-{\mathbf{d}}_{\mathrm{q}}(\mathbf{q},{\mathbf{q}}_0)\right)+\left(1-s\right){\mathbf{d}}_{\mathrm{q}}({\mathbf{q}}_{\mathrm{g}},\mathbf{q})\end{array}\right]-\left[\begin{array}{cc}{\mathbf{D}}_{\mathrm{p}}& \\ & {\mathbf{D}}_{\mathrm{q}}\end{array}\right]\left[\begin{array}{c}{\mathbf{v}}_{\mathrm{d}}\\ {\mathit{\omega }}_{\mathrm{d}}\end{array}\right],$$

(18)

with

$$\tau \left[\begin{array}{c}\dot{\mathbf{x}}\\ \dot{\mathbf{q}}\end{array}\right]=\left[\begin{array}{c}\mathbf{v}\\ \mathit{\eta }*\mathbf{q}/2\end{array}\right],\tau \left[\begin{array}{c}{\mathbf{a}}_{\mathrm{d}}\\ {\mathit{\alpha }}_{\mathrm{d}}\end{array}\right]=\left[\begin{array}{c}{\dot{\mathbf{v}}}_{\mathrm{d}}\\ {\dot{\mathit{\omega }}}_{\mathrm{d}}\end{array}\right]$$

(19)

where ${\mathbf{K}}_{\mathrm{p}},{\mathbf{K}}_{\mathrm{q}},{\mathbf{D}}_{\mathrm{p}},{\mathbf{D}}_{\mathrm{q}}$ are constant coefficients. $\tau$ is the time scaling factor. $\mathbf{x}$ denotes the end-effector position, ${\mathbf{x}}_0$ the initialing position, and ${\mathbf{x}}_{\mathrm{g}}$ the final position. $\mathbf{q}$ denotes the quaternion of the end effector’s orientation, ${\mathbf{q}}_0$ the initialing quaternion, and ${\mathbf{q}}_{\mathrm{g}}$ the final quaternion. $\mathbf{v}$ denotes the velocity of the end effector and $\mathit{\omega }$ the angular velocity. ${\mathbf{a}}_{\mathrm{d}}$ denotes the acceleration of the end effector and ${\mathit{\alpha }}_{\mathrm{d}}$ the angular acceleration. $\mathit{\eta }$ is a quaternion with a scalar part of 0 and an imaginary part of $\mathit{\omega }$. ${\mathbf{d}}_{\mathrm{q}}:{\mathbb{R}}^4\mapsto \mathrm{so}\left(3\right)$ is the logarithmic mapping of the quaternion deviation to angular velocity space, as defined in [28]. $\mathbf{g}$ denotes the radial basis vector, and ${\mathbf{w}}_{\mathrm{p}}$ and ${\mathbf{w}}_{\mathrm{q}}$ are the parameter matrices of $N\times 3$. The $n\mathrm{th}$ element of $\mathbf{g}$ is defined as follows

$${\left[\mathbf{g}\right]}_n=\frac{{\varphi }_n\left(s\right)}{{\sum }_{n=1}^N{\varphi }_n\left(s\right)},$$

(20)

with,

$${\varphi }_n\left(s\right)=exp(-h_n{(s-c_n)}^2),$$

(21)

where ${\varphi }_n\left(s\right)$ denotes the $n\mathrm{th}$ Gaussian radial basis function, $h_n$ the width of this basis, and $c_n$ the center of this basis.

The Equation (18) defines a virtual critical damping system with attraction points ${[{\mathbf{x}}_{\mathrm{g}},{\mathbf{q}}_{\mathrm{g}}]}^{\mathrm{T}}$ and ${[{\mathbf{w}}_{\mathrm{p}}^{\mathrm{T}}\mathbf{g},{\mathbf{w}}_{\mathrm{q}}^{\mathrm{T}}\mathbf{g}]}^{\mathrm{T}}$. ${[{\mathbf{x}}_{\mathrm{g}},{\mathbf{q}}_{\mathrm{g}}]}^{\mathrm{T}}$ converges the system towards the final position and ${[{\mathbf{w}}_{\mathrm{p}}^{\mathrm{T}}\mathbf{g},{\mathbf{w}}_{\mathrm{q}}^{\mathrm{T}}\mathbf{g}]}^{\mathrm{T}}$ adjusts the shape of the motion trajectory. The importance of ${[{\mathbf{x}}_{\mathrm{g}},{\mathbf{q}}_{\mathrm{g}}]}^{\mathrm{T}}$ and ${[{\mathbf{w}}_{\mathrm{p}}^{\mathrm{T}}\mathbf{g},{\mathbf{w}}_{\mathrm{q}}^{\mathrm{T}}\mathbf{g}]}^{\mathrm{T}}$ to the system is weighted by $1-s$ and s. When t is small, the shape of the generated trajectory is more important, and as s trends toward 0, convergence to the final position ${[{\mathbf{x}}_{\mathrm{g}},{\mathbf{q}}_{\mathrm{g}}]}^{\mathrm{T}}$ becomes more important.

(3) Learning Parameters: ${\mathbf{Y}}_{\mathrm{p}}$ and ${\mathbf{Y}}_{\mathrm{q}}$ are the targets derived from the demonstration trajectory time series,

$$\begin{array}{cc}\hfill & {\mathbf{Y}}_{\mathrm{p}}=[{\mathbf{y}}_{\mathrm{p}}\left(t_1\right),{\mathbf{y}}_{\mathrm{p}}\left(t_2\right),\dots ,{\mathbf{y}}_{\mathrm{p}}\left(t_M\right)],\hfill \\ \hfill & {\mathbf{Y}}_{\mathrm{q}}=[{\mathbf{y}}_{\mathrm{q}}\left(t_1\right),{\mathbf{y}}_{\mathrm{q}}\left(t_2\right),\dots ,{\mathbf{y}}_{\mathrm{q}}\left(t_M\right)],\hfill \end{array}$$

(22)

with

$$\begin{array}{cc}\hfill {\mathbf{y}}_{\mathrm{p}}\left(t\right)& =-s[{\mathbf{x}}_{\mathrm{ref}}\left(t\right)-{\mathbf{x}}_{\mathrm{ref}}\left(t_1\right)]+(1-s)[{\mathbf{x}}_{\mathrm{ref}}\left(t_M\right)-{\mathbf{x}}_{\mathrm{ref}}\left(t\right)-\frac{{\mathbf{D}}_{\mathrm{p}}{\mathbf{v}}_{\mathrm{ref}}\left(t\right)}{{\mathbf{K}}_{\mathrm{p}}}-\frac{{\mathbf{a}}_{\mathrm{ref}}\left(t\right)}{{\mathbf{K}}_{\mathrm{p}}},\hfill \\ \hfill {\mathbf{y}}_{\mathrm{q}}\left(t\right)& =-s\mathbf{d}({\mathbf{q}}_{\mathrm{ref}}\left(t\right),{\mathbf{q}}_{\mathrm{ref}}\left(t_1\right))+\left(1-s\right)\mathbf{d}({\mathbf{q}}_{\mathrm{ref}}\left(t_M\right),{\mathbf{q}}_{\mathrm{ref}}\left(t\right))-\frac{{\mathbf{D}}_{\mathrm{q}}{\mathit{\omega }}_{\mathrm{ref}}\left(t\right)}{{\mathbf{K}}_{\mathrm{q}}}-\frac{{\mathit{\alpha }}_{\mathrm{ref}}\left(t\right)}{{\mathbf{K}}_{\mathrm{q}}}.\hfill \end{array}$$

(23)

The reference trajectory and feedforward force jointly describe the cutting behavior. The parameter matrices of the nonlinear terms ${\mathbf{w}}_{\mathrm{p}}$ and ${\mathbf{w}}_{\mathrm{q}}$ are the representations of the cutting behavior in the parameter space.

${\mathbf{w}}_{\mathrm{p}}$ and ${\mathbf{w}}_{\mathrm{q}}$ are learned from the demonstration data by solving the optimization problems (24),

$$\begin{array}{cc}\hfill min& \sum _{i=1}^M\parallel {\mathbf{w}}_{\mathrm{p}}^{\mathrm{T}}\mathbf{g}\left(s_i\right)s_i-{\mathbf{y}}_{\mathrm{p}}\left(t_i\right)\parallel ,\hfill \\ \hfill min& \sum _{i=1}^M\parallel {\mathbf{w}}_{\mathrm{q}}^{\mathrm{T}}\mathbf{g}\left(s_i\right)s_i-{\mathbf{y}}_{\mathrm{q}}\left(t_i\right)\parallel ,\hfill \end{array}$$

(24)

where, at each time step, i, $s_i$ represents the activation signal, and $\mathbf{g}\left(s_i\right)$ is the corresponding radial basis vector. The parameter matrices to be solved are denoted as ${\mathbf{w}}_{\mathrm{p}}$ and ${\mathbf{w}}_{\mathrm{q}}$. We utilized the locally weighted linear regression (LWR) method to solve Equation (24) in this study.

2.4.2. Robot Joints Admittance Control

In this study, we implement the proposed algorithm on an Aubo-i7 robot arm (Shanghai, China). The arm has six joints. The controller interacts with the low-level embedded motion controller of the robot to set and read the joints’ position. The interaction is through robot operating system (ROS) software architecture and TCP-CAN transparent transmission at 50 Hz. Regarding the robot body as rigid, we design the following admittance control law (25) to the joints’ motion for safe and effective cutting.

$${\mathbf{K}}_{\mathrm{j}}({\mathit{\theta }}_{\mathrm{m}}-{\mathit{\theta }}_{\mathrm{r}})-{\mathbf{D}}_{\mathrm{j}}({\dot{\mathit{\theta }}}_m-{\dot{\mathit{\theta }}}_r)=\Delta \mathit{\tau },$$

(25)

where ${\mathbf{k}}_{\mathrm{j}}$ is the joint stiffness coefficient matrix and ${\mathbf{d}}_{\mathrm{j}}$ the joint damping coefficient matrix. ${\mathit{\theta }}_{\mathrm{m}}$ is the command positions of the joints. $\delta \mathit{\tau }$ is the joints’ torque offsets. We conduct the following transformation (26) to transform the cutting force offsets in the coordinate system of the end effector into the joints’ torque offsets

$$\Delta \mathit{\tau }={\mathbf{J}}^{\mathrm{T}}({\mathbf{F}}_{\mathrm{fb}}-{\mathbf{F}}_{\mathrm{ff}}),$$

(26)

where, ${\mathbf{F}}_{\mathrm{ff}}$ is the feedforward force, ${\mathbf{F}}_{\mathrm{fb}}$ is the measured feedback force. ${\mathbf{K}}_{\mathrm{j}}$ and ${\mathbf{D}}_{\mathrm{j}}$ are calculated by Equations (27) and (28)

$$\begin{array}{cc}\hfill {\mathbf{K}}_{\mathrm{j}}& ={\mathbf{J}}^{\mathrm{T}}\mathbf{\Gamma }\mathbf{KJ},\hfill \end{array}$$

(27)

$$\begin{array}{cc}\hfill {\mathbf{D}}_{\mathrm{j}}& ={\mathbf{J}}^{\mathrm{T}}{\mathbf{\Gamma }}^{\frac{1}{2}}\mathbf{D}\mathbf{J},\hfill \end{array}$$

(28)

where $\mathbf{D}=2\sqrt{\mathbf{K}\mathbf{M}}$, $\mathbf{M}=\left[\begin{array}{ccc}{m}_1& \cdots & 0\\ ⋮& \ddots & ⋮\\ 0& \cdots & m_6\end{array}\right]$, $\mathbf{K}=\left[\begin{array}{ccc}{k}_1& \cdots & 0\\ ⋮& \ddots & ⋮\\ 0& \cdots & k_6\end{array}\right]$, and $\mathbf{\Gamma }=\left[\begin{array}{ccc}{\gamma }_1& \cdots & 0\\ ⋮& \ddots & ⋮\\ 0& \cdots & {\gamma }_6\end{array}\right]$. $\mathbf{k}$ is the stiffness coefficient matrix on the end effector coordinates, and $\mathbf{d}$ is the damping coefficient matrix on the end effector coordinates. $\mathbf{\Gamma }$ is the stiffness adjustment factor matrix, and its diagonal elements ${\gamma }_1,\dots ,{\gamma }_6$ are the stiffness adjustment factors in the direction of the corresponding axis of the end effector. The initial values of ${\gamma }_1,\dots ,{\gamma }_6$ are all 1. The stiffness and damping of the admittance control in different directions can be set separately by changing ${\gamma }_1,\dots ,{\gamma }_6$.

The admittance control converts the end effector force offsets into the angular velocity of the robot’s joints. Based on the Equations (25)–(28), we propose the following update law (29) and (30) of the velocities and angles of the joints

$$\begin{array}{cc}\hfill {\dot{\mathit{\theta }}}_{\mathrm{m}}& \leftarrow {\dot{\mathit{\theta }}}_{\mathrm{r}}+{\left({\mathbf{J}}^{\mathrm{T}}\mathbf{D}\mathbf{J}\right)}^{-1}\left[{\mathbf{J}}^{\mathrm{T}}({\mathbf{F}}_{\mathrm{fb}}-{\mathbf{F}}_{\mathrm{ff}})-{\mathbf{J}}^{\mathrm{T}}\mathbf{KJ}({\mathit{\theta }}_{\mathrm{m}}-{\mathit{\theta }}_{\mathrm{r}})\right],\hfill \end{array}$$

(29)

$$\begin{array}{cc}\hfill {\mathit{\theta }}_{\mathrm{m}}& \leftarrow {\mathit{\theta }}_{\mathrm{m}}+{\dot{\mathit{\theta }}}_{\mathrm{m}}\mathrm{d}t,\hfill \end{array}$$

(30)

After that, the joints’ angles are sent to the low-level embedded motion controller of the robot for execution. This paper does not include a discussion of the low-level embedded motion controller.

2.4.3. Inverse Velocity Motion Control

The robot updates the reference pose and feedforward force each time step during cutting. First, it updates ${\dot{\mathbf{v}}}_{\mathrm{d}}$ and ${\dot{\mathit{\omega }}}_{\mathrm{d}}$ using Formula (18). Moreover, ${\mathbf{v}}_d$ and ${\mathit{\omega }}_d$ are updated via numerical integrals (31)

$$\left[\begin{array}{c}{\mathbf{v}}_{\mathrm{d}}\\ {\mathit{\omega }}_{\mathrm{d}}\end{array}\right]\leftarrow \left[\begin{array}{c}{\mathbf{v}}_{\mathrm{d}}\\ {\mathit{\omega }}_{\mathrm{d}}\end{array}\right]+\left[\begin{array}{c}{\dot{\mathbf{v}}}_{\mathrm{d}}\\ {\dot{\mathit{\omega }}}_{\mathrm{d}}\end{array}\right]\Delta t,$$

(31)

Then, the reference position ${\mathbf{x}}_d$ and the reference pose ${\mathbf{q}}_d$ are updated

$$\begin{array}{cc}\hfill {\mathbf{x}}_d& \leftarrow \mathbf{x}+\frac{{\mathbf{v}}_d\Delta t}{\tau }\hfill \end{array}$$

(32)

$$\begin{array}{cc}\hfill {\mathbf{q}}_d& \leftarrow {exp}_{\mathrm{q}}\left(\frac{{\mathit{\omega }}_d\Delta t}{2\tau }\right)*\mathbf{q}.\hfill \end{array}$$

(33)

To achieve continuous tracking of the end effector’s motion and avoid joint angle and velocity limits, we need to solve for the inverse kinematics [29,30,31] and obtain the command joint angles ${\mathit{\theta }}_{\mathrm{r}}$, as described in Equation (25). In detail, we use the instantaneous kinetic energy of the robot joints as the cost function. The constraints we take into account include the first-order forward kinematics of the robot and the physical limitations of the joint angles and velocities. The joint velocities are then obtained by solving the following optimization problem

$$\begin{array}{cc}\hfill \underset{{\dot{\mathit{\theta }}}_{\mathrm{r}}}{min}& {\dot{\mathit{\theta }}}_{\mathrm{r}}^{\mathrm{T}}\mathbf{M}{\dot{\mathit{\theta }}}_{\mathrm{r}},\hfill \end{array}$$

(34)

$$\begin{array}{cc}\hfill \mathrm{s}.\mathrm{t}.& \mathbf{J}{\dot{\mathit{\theta }}}_{\mathrm{r}}={[{\mathbf{v}}_{\mathrm{r}},{\mathit{\omega }}_{\mathrm{r}}]}^{\mathrm{T}},\hfill \end{array}$$

(35)

$$\begin{array}{cc}\hfill & {\mathit{\theta }}_{\mathrm{r}}^{min}\le {\mathit{\theta }}_{\mathrm{r}}\le {\mathit{\theta }}_{\mathrm{r}}^{max},\hfill \end{array}$$

(36)

$$\begin{array}{cc}\hfill & {\dot{\mathit{\theta }}}_{\mathrm{r}}^{min}\le {\dot{\mathit{\theta }}}_{\mathrm{r}}\le {\dot{\mathit{\theta }}}_{\mathrm{r}}^{max}.\hfill \end{array}$$

(37)

where $\mathbf{M}$ is the inertial matrix. Equation (34) is the instantaneous kinetic energy of the joint to be optimized. Equation (35) is the forward kinematic constraint. Equation (36) is the constraint on the physical position of the joint angle. Equation (37) is the physical velocity constraint of the joint. The solution to the above optimization problem is

$${\dot{\mathit{\theta }}}_{opt}={\mathbf{J}}^{†}{[{\mathbf{v}}_{\mathrm{d}},{\mathit{\omega }}_{\mathrm{d}}]}^{\mathrm{T}}.$$

(38)

In Equation (38), the pseudo-inverse ${\mathbf{J}}^{†}$ of the joint Jacobi $\mathbf{J}$ is calculated as

$${\mathbf{J}}^{†}={\mathbf{M}}^{-1}{\mathbf{J}}^{\mathrm{T}}\left(\mathbf{J}{\mathbf{M}}^{-1}{\mathbf{J}}^{\mathrm{T}}\right).$$

(39)

In this study, we have introduced compensation in Formula (38) to eliminate the cumulative error of numerical calculation. It ensures that the robot tracks the motions generated by the modulated DMPs. The final calculation Formula (40) of ${\dot{\mathit{\theta }}}_r$ is as follows:

$${\dot{\mathit{\theta }}}_r={\dot{\mathit{\theta }}}_{opt}+{\mathbf{J}}^{†}{log}_6\left({\mathbf{X}}_r^{-1}{\mathbf{X}}_d\right),$$

(40)

where ${\mathbf{X}}_r$ is the homogeneous matrix of the end effector pose at ${\mathit{\theta }}_r$ derived from the forward kinematics of the robot. ${\mathbf{X}}_d$ is the homogeneous matrix of the expected end effector pose, which corresponds to ${\mathbf{x}}_d,{\mathbf{q}}_d$. ${log}_6:\mathbf{SE}\left(3\right)\mapsto {\mathbb{R}}^6$ is the logarithmic mapping of the robot’s homogeneous matrix to the twist in motion.

3. Results and Discussion

3.1. Human–Robot Skill Transfer

We train the proposed method on the datasets mentioned in Section 2.2 to demonstrate its effectiveness and applicability in learning and generating soft tissue cutting behaviors.

It can be observed from Figure 4 that the PCA eigenvector heat maps vary remarkably in scenarios 1, 2, 3, and 4. This suggests that the latent space variables in Equation (6) are able to effectively identify the cutting behaviors of four scenarios, although we use PCA to calculate one-dimensional latent space variables for simplicity.

Table 1. Results of latent variable encoding and generation.

Scenario	Samples	w1	w2	w3	w4	w5	w6	PEavg (m)	OEavg (rad)	Ratio (%)
1	4	0.3864	−0.01621	−0.3109	−0.4989	0.4202	−0.5730	0.056	0.032	55.50
2	2	0.4390	−0.06339	0.1329	−0.8468	0.2571	−0.04962	0.038	0.025	80.96
3	2	−0.6637	−0.1099	0.6496	−0.04189	−0.3370	−0.1010	0.030	0.040	48.80
4	30	−0.5927	−0.2432	−0.4050	−0.4467	−0.2707	0.3908	0.047	0.038	74.71

$w_1$, $w_2$, $w_3$, $w_4$, $w_5$, and $w_6$ denote the elements of the PCA eigenvector.

shows the explanation rates of GMM, axis weights of PCA, and pose reconstruction errors of GMR in four scenarios. We find in Table 1 that Scenario 1 has the smallest difference in axis weights, indicating a more complex contact situation. In this scenario, there are both contacts between the side of the cutting blade and the ischial surface and contact between the cutting blade and the ischial surface. Moreover, it can be observed from Table 1 that $w_2$ is small among the four scenarios. This indicates that when the side of the cutting knife is in contact with the rigid medium, the large elastic deformation may reduce the measured value of the cutting force in this direction. From the perspective of the interpretation rate, on the one hand, the interpretation rates of Scenarios 4 and 2 are 74.71% and 80.96%, indicating that in these two scenarios, the GMM-GMR module in this work can effectively represent the cutting behavior in the latent space. On the other hand, the explanation rates of scenarios 1 and 3 are nearly 50.00%, indicating that in these two scenarios, the effectiveness of latent space variables in representing cutting behavior can be further improved.

In addition, in Scenarios 2, 3, and 4, the position reconstruction error is less than 0.05 m, the position reconstruction error of Scenario 3 is close to 0.05 m, and in the four scenarios, the attitude reconstruction error is less than 0.05 rad. Considering the jitter and operating accuracy during human demonstrations, this suggests that the latent space GMM-GMR module can effectively generate the target cutting behaviors for the robot in all the four scenarios.

3.2. Foam Cutting Test

We carry out a foam-cutting experiment to verify the reliability and effectiveness of the proposed framework in adapting to the position changes of the rigid medium under soft tissue. The robot uses the learned parameters and fixed stiffness coefficient to perform one cutting each at the loading platform heights of +30 mm, +20 mm, +10 mm, 0 mm, −10 mm, −20 mm, and −30 mm. The to-be-cut soft tissue (a foam block) is pasted on the surface of the rigid medium, which is a plastic pad used for food processing. We fix the plastic pad on the loading robot platform, and the loading robot platform can adjust the height of the plastic pad by itself. The size of the foam block used is 60 × 50 × 50 mm, and the height of the plastic backing plate at the beginning of the test is 0 mm. As is shown in the Figure 5, the cutting process can be divided into four phases, i.e., entering, cutting, releasing, and returning.

As shown in Figure 6, in each test, the robot successfully separates the foam. The first phase is 0–5.0 s, the second phase is 5.0–10.5 s, the third phase is 10.5–15.0 s, and the fourth phase is 15.0–20.0 s. In the first phase, the main contact direction is the Z-axis. The direction of $F_{\mathrm{x}}$ is negative, the direction of $F_{\mathrm{y}}$ is positive, and the direction of $F_{\mathrm{z}}$ is negative. On the one hand, the changes in $F_{\mathrm{x}}$ and $F_{\mathrm{y}}$ are small, falling from 0 N to approximately −5 N and rising from 0 N to approximately 5 N, respectively. Both $F_{\mathrm{x}}$ and $F_{\mathrm{y}}$ are within the safe range of the force sensor and robot payload. On the other hand, $F_{\mathrm{z}}$ decreases rapidly. The second phase is the cutting phase. $F_{\mathrm{y}}$ rises from approximately 5 N to 25 N and maintains the great pressure within the safety range of the force sensor and robot load. In this phase, $F_{\mathrm{y}}$ is mainly used to apply pressure on the plastic plate to ensure effective cutting of the thin layer in the boundary between the foam and the plastic plate and to compensate for the resistance of the plastic backing plate to the cutting tool. $F_{\mathrm{x}}$ further dropped from −5 N to −10 N. $T_{\mathrm{y}}$ and $F_{\mathrm{x}}$ have similar trends. Both $F_{\mathrm{x}}$ and $T_{\mathrm{y}}$ are close to the data during human demonstration. $F_{\mathrm{x}}$ and $T_{\mathrm{y}}$ are mainly derived from the demonstration. When demonstrating, the plane of the blade and the plastic plate are not perpendicular. The third stage is the release stage. The cutting force of each axis gradually returns to 0 N. From the perspective of motion, there is only a minor deviation between the test trajectories and the demonstration trajectories. This shows that after the cutting blade leaves the plastic plate and foam, the IVAC algorithm can effectively adjust the motion trajectories and reduce the offset between the target trajectory and the actual trajectory.

Table 2. Foam cutting forces during robot cutting at different platform heights.

		+30 mm	+20 mm	+10 mm	+00 mm	−10 mm	−20 mm	−30 mm	Human
	maximum	2.459	1.646	0.285	0.233	0.245	0.248	0.317	0.682
$F_{\mathrm{x}}$	minimum	−9.376	−9.301	−9.054	−8.683	−8.583	−8.548	−8.235	−10.271
	average	−2.560	−2.582	−2.579	−2.383	−2.286	−2.227	−2.158	−3.177
	maximum	34.400	33.778	32.592	31.604	30.917	30.180	29.119	36.452
$F_{\mathrm{y}}$	minimum	−0.063	−0.155	0.001	−0.031	−0.031	−0.133	−0.237	−0.189
	average	12.463	11.747	11.136	10.558	10.200	9.749	9.202	12.985
	maximum	20.068	20.783	13.843	12.388	13.660	15.097	13.447	13.313
$F_{\mathrm{z}}$	minimum	−21.865	−19.712	−22.068	−19.825	−16.831	−13.500	−12.252	−11.229
	average	2.289	2.698	1.468	0.923	1.467	2.161	2.039	2.084
	maximum	0.028	0.004	0.004	0.003	0.003	0.004	0.008	0.016
$T_{\mathrm{x}}$	minimum	−9.282	−9.112	−8.871	−8.682	−8.548	−8.392	−8.178	−8.478
	average	−3.237	−3.105	−2.940	−2.827	−2.741	−2.658	−2.540	−2.747
	maximum	0.727	0.485	0.046	0.044	0.044	0.048	0.061	0.041
$T_{\mathrm{y}}$	minimum	−2.371	−2.335	−2.315	−2.280	−2.259	−2.249	−2.215	−2.241
	average	−0.653	−0.643	−0.635	−0.619	−0.607	−0.598	−0.571	−0.657
	maximum	0.003	0.003	0.002	0.006	0.004	0.002	0.003	0.001
$T_{\mathrm{z}}$	minimum	−0.227	−0.227	−0.204	−0.217	−0.213	−0.204	−0.192	−0.202
	average	−0.047	−0.040	−0.036	−0.036	−0.038	−0.037	−0.033	−0.041

presents the maximum, minimum, and average forces on each axis of the force sensor during robot cutting at different platform heights, as well as the reference values measured during human teaching. When the platform height is lower than −10 mm, its amplitude is close to the human demonstration. As the platform height increases, The amplitude of $F_{\mathrm{z}}$ is also larger. However, in each test, $F_{\mathrm{z}}$ is still within the safe range of the force sensor and robot load.

The above results suggest that the proposed IVAC algorithm can adapt to different platform heights quickly by changing its motion trajectory and has a wide adjustment range. The proposed framework can ensure that the cutting force is within a safe range, which is suitable for cutting when the position of the rigid medium under soft tissue is uncertain. It can also be found from the results that the cutting force errors of the same axis are quite different at different execution phases. It suggests that the optimal stiffness differs from phases and axes during execution. Follow-up research can consider the optimization of the stiffness coefficients of IVAC through reinforcement learning or variable admittance control.

3.3. Sheep Hindquarters Separation Test

This study carries out a sheep hindquarter separation test to evaluate the effectiveness of the proposed framework in the meat cutting scenario. The separation processes of the left hind leg of the hindquarter are selected as the test scenarios. In the test, the robot needed to cut the leg muscles along the hip bone (ilium and ischium) of the hindquarter to separate the left hind-leg from the hip bone. To complete the separation, the robot uses the proposed IVAC algorithm and the learned skill to perform cutting in Scenarios 1–3 sequentially. We prepare two sheep hindquarters for testing.

The separation effects are shown in Figure 7; the robot successfully separates the hip bone and hind-leg of the sheep carcass. The hind-leg of the divided sheep carcass is relatively complete, and there is only a little meat left on the hip bone. Moreover, no bone cracks appeared during the separation processes, and there is no collision accident between the cutting tool and the bone. This indicates that the cutting by the robot is compliant and accurate.

Table 3. Performance of sheep carcass hindquarter cutting test from 2 samples.

Criete	Sample 1			Sample 2
	Scenario 1	Senario 2	Scenario 3	Scenario 1	Scenario 2	Scenario 3
Weight [kg]	5.18			6.30
$RR$ [%]	5.6			4.8
Maximal thickness [mm]	3.1			3.8
${\left|F_{\mathrm{x}}\right|}_{\max }$[N]	6.461	3.056	4.072	7.655	7.441	4.787
${\left|F_{\mathrm{y}}\right|}_{\max }$[N]	30.719	24.780	15.972	31.981	27.404	20.432
${\left|F_{\mathrm{z}}\right|}_{\max }$[N]	17.968	21.855	4.185	30.686	22.858	5.004
${\left|T_{\mathrm{x}}\right|}_{\max }$[N]	5.708	4.971	3.803	5.726	5.2870	4.808

shows the meat residual rate (RR), maximum residual meat thickness, and maximum absolute value of cutting force. The maximum residual meat thickness is used to characterize the adaptive performance of IVAC. The smaller the maximum residual meat thickness, the better the adaptive performance of IVAC when there are differences in shape and size between the test sample and the teaching sample. The maximum absolute value of the cutting force reflects the cutting resistance under the premise of completing the cutting. The smaller the value, the higher the cutting efficiency. Moreover, RR is defined as the ratio of the weight of the residual hind leg meat on the hip bone of the sheep carcass to the total weight of the separated hind-leg. RR evaluates the loss of high-value leg meat during the cutting process. The smaller the value, the less the loss. RR is calculated with Equation (41)

$$RR=\frac{2.0\times m_{\mathrm{res}}}{m_{\mathrm{res}}+m_{\mathrm{leg}}}\times 100\%,$$

(41)

where $RR$ denotes the residual rate, $m_{\mathrm{res}}$ denotes the weight of the residual hind-leg meat on the hip bone, and $m_{\mathrm{leg}}$ denotes the weight of the separated hind-leg.

According to Table 3, the maximum residual meat thickness of sample 1 is 3.1 mm and sample 2 is 3.8 mm, both of which are less than 5 mm. The segmentation residual rates of samples 1 and 2 are 5.6% and 4.8%. The results indicate that the proposed method can segment sheep carcasses well and has good generalization and effectiveness. In addition, the maximum residual meat thickness and segmentation residual rate obtained in the experiment are small, indicating that the segmentation accuracy is high enough to segment the hind-legs of sheep carcasses. The last four rows of Table 3 are the absolute maximum value of the cutting force in the directions of $F_{\mathrm{x}}$, $F_{\mathrm{y}}$, $F_{\mathrm{z}}$, and $T_{\mathrm{x}}$. It can be observed from Table 3 that the largest absolute values of the cutting force in each axis of samples 1 and 2 under the three scenarios are all within the safe range of the robot and the force sensor. This suggests that the proposed method can not only ensure that the cutting depth and residual rate meet the process requirements but also avoid excessive cutting force through adaptation. The proposed framework can realize the human–robot skill transfer in the sheep hindquarters separation test.

Figure 8 shows the cutting forces and torque curves of cutting sample 1 to illustrate the similarity between robot cutting and the human demonstration. The first column of Figure 8 shows the results of cutting scenario 1. (a) The robotic cutting force tendency of each axis is consistent with the reference curve measured in the human demonstration, indicating that the robot can imitate the human cutting behavior well. (b) The robot cutting force measured is close to the humans demonstration. (c) The maximum robot cutting force in the $F_{\mathrm{z}}$ direction (200 steps) is greater than the human cutting force. However, it is still within the safe range of the robot and force sensor. This suggests that the proposed method can perform cutting efficiently and avoid excessive cutting force.

The values in the second column of Figure 8 are curves in scenario 2. Meanwhile, the tendency and amplitude of the robot’s cutting force $F_{\mathrm{y}},F_{\mathrm{z}},T_{\mathrm{x}}$ are very close to that of human cutting. As is shown in the third column of Figure 8, the tendency and amplitude of the robot cutting force are near human, and the vibration of robot cutting is greater in scenario 3. Additionally, Figure 8 shows that the total separation time of the hind-leg is nearly 30 s. This suggests that the proposed method has practical efficiency compared with manual cutting.

4. Conclusions

The meat industry requires automation due to labor shortages and inefficiencies. Traditional 2D methods are limited when cutting complex 3D anatomy structures, such as sheep hindquarters. Compliant cutting in a 3D space can precisely separate the meat tissue from the bones but lacks flexibility. A multi-demonstration HRST framework is proposed to address these challenges, reduce cutting force, and minimize meat loss. This method potentially benefits flexible and generalized robot skill learning in interaction tasks, which provides an alternative approach for harvesting high-value cuts and increasing yields in meat industries. Multi-demonstrations of cutting are encoded through the latent space GMM-GMR module to generate the target cutting behaviors for robot, which enables the generalization to meat parts and sizes. The IVAC-DMP control provides the flexibility of modulating robot joint motion during cutting and avoids excessive cutting forces.

The experimental results show that the latent space GMM-GMR can effectively extract cutting skills from multi-demonstrations. The foam-cutting test results show that the IVAC module can adapt to the position changes of the rigid medium under soft tissue reliably and efficiently. The sheep hindquarters separation test result shows that the maximum residual meat thickness of sample 1 and 2 are 3.1 mm and 3.8 mm. The residual rates of samples 1 and 2 are 5.6% and 4.8%. Both meet the requirements of the sheep hindquarter separation.

There are several limitations to the proposed framework. First, the current PCA representation of cutting behavior may not be completely accurate. In the future, it may be beneficial to use models with stronger representation capabilities, e.g., variational autoencoder to improve the accuracy of the latent space representation. Second, the demonstrating scenarios used in this study are limited to cutting sheep hindquarters and foam blocks. To extract more generalizable semantic cutting skills, future studies could apply the framework to broader tissue-cutting tasks by training with larger-scale datasets from various meat pieces over longer epochs. Finally, the IVAC module utilized fixed admittance parameters. To improve interaction flexibility, future studies could apply variable AC or IC to provide optimal stiffness at different stages of the cutting process and in different cutting directions.