A Trajectory Optimisation-Based Incremental Learning Strategy for Learning from Demonstration

Abstract

The insufficient generalisation capability of the conventional learning from demonstration (LfD) model necessitates redemonstrations. In addition, retraining the model can overwrite existing knowledge, making it impossible to perform previously acquired skills in new application scenarios. These are not economical and efficient. To address the issues, in this study, a broad learning system (BLS) and probabilistic roadmap (PRM) are integrated with dynamic movement primitive (DMP)-based LfD. Three key innovations are proposed in this paper: (1) segmentation and extended demonstration: a 1D-based topology trajectory segmentation algorithm (1D-SEG) is designed to divide the original demonstration into several segments. Following the segmentation, a Gaussian probabilistic roadmap (G-PRM) is proposed to generate an extended demonstration that retains the geometric features of the original demonstration. (2) DMP modelling and incremental learning updating: BLS-based incremental learning for DMP (Bi-DMP) is performed based on the constructed DMP and extended demonstration. With this incremental learning approach, the DMP is capable of self-updating in response to task demands, preserving previously acquired skills and updating them without training from scratch. (3) Electric vehicle (EV) battery disassembly case study: this study developed a solution suitable for EV battery disassembly and established a decommissioned battery disassembly experimental platform. Unscrewing nuts and battery cell removal are selected to verify the effectiveness of the proposed algorithms based on the battery disassembly experimental platform. In this study, the effectiveness of the algorithms designed in this paper is measured by the success rate and error of the task execution. In the task of unscrewing nuts, the success rate of the classical DMP is 57.14% and the maximum error is 2.760 mm. After the optimisation of 1D-SEG, G-PRM, and Bi-DMP, the success rate of the task is increased to 100% and the maximum error is reduced to 1.477 mm.

1. Introduction

Industrial robots have been designed to perform repetitive tasks, so they have been extensively used in the mass production of automotive, aerospace, and marine industries. As the demand for customised production has dramatically increased in recent years, collaborative robots (cobots) have been actively developed to carry out customised tasks within unstructured environments [1]. Learning from demonstration (LfD) is an effective approach to facilitating cobots to fulfil the above purpose. LfD can enable cobots to efficiently realise customised operations based on human demonstrations without the need for detailed motion programming or reprogramming [2,3]. The process of LfD is illustrated in Figure 1, where a cobot performs a pick-and-place task. To create a demonstration, a human operator drags the end-effector of the cobot to pick up an object from a start point (the yellow box in Figure 1) and place it at a target point (the marker in Figure 1). The demonstration (trajectory) is shown as the red curve in Figure 1. The demonstration is recorded and used to train an LfD-based learning model. To accomplish a similar task where the marker is moved to a new position, the trained LfD-based learning model can adaptively generate a new trajectory (generated trajectory). This is illustrated as the blue dashed line in Figure 1.

Dynamic movement primitive (DMP) is an effective learning model for implementing LfD [4]. Unlike Gaussian mixture modelling (GMM), the DMP model only needs to learn from one demonstration, thereby minimising the difficulty of creating multiple demonstrations [5,6,7]. The mathematical representation of the DMP model includes a second-order spring damping system and a force item [8]. The former ensures that the generated trajectory created using the DMP model can converge to the target without divergence. The latter, which is constructed based on a series of radial basis function (RBF) kernels, is used to control the shape of the generated trajectory and manage the convergence process. The advantage of the DMP model is its robustness to perturbations [9]. However, the model exhibits some limitations as well. For instance, when the position of the marker illustrated in Figure 1 is moved far away from the previous demonstration area, the generated trajectory via DMP cannot work properly [10]. That is, the difference between the positions of the generated trajectory and the marker should not exceed δ, which is a threshold representing the maximum error allowable for the successful execution of the pick-and-place task. The situation is illustrated by the green dashed curve in Figure 1. To address the issue, the following considerations are given:

• A straightforward idea is to create a new demonstration for the new situation. However, redemonstration is time consuming as the experimental environment needs to be set up again [3]. To mitigate the problem, it will be beneficial to reuse part of the originally created demonstration to adaptively generate a new trajectory (namely, an extended demonstration) for the new target;
• DMP is a one-shot learning model, which means that when the extended demonstration is learned, the knowledge gained from the previous demonstration is forgotten. Therefore, it is imperative to develop an incremental learning mechanism for the DMP model to improve its generalisation capability for various situations.

Various approaches for trajectory generation were developed to complement demonstration data, such as task-parameterised Gaussian mixture models (TP-GMMs) [3], adversarial generative models [11], B-spline-based models [12], and PRM-based models [13]. These approaches can generate new trajectories flexibly. However, the lack of imitation and utilisation of prior knowledge (i.e., existing demonstrations) poses safety risks in LfD applications. In addition, some researchers attempted to improve the performance of DMP [14]. For instance, neural networks and reinforcement learning [15,16,17] were used to enhance the learning performance of the force term in the DMP model. The development of the incremental learning mechanisms for the DMP model received attention as well [18,19,20]. However, the previous research has problems of high complexity and computational difficulty. A more detailed analysis is provided in the following Section 2.

Based on the above considerations, in this study, an improved DMP model with extended demonstration and incremental learning capabilities is designed. The capabilities are implemented based on a 1D-based topology trajectory segmentation algorithm (1D-SEG), a Gaussian probabilistic roadmap (G-PRM), and a broad learning system (BLS). By constructing an appropriate sampling strategy and loss function, G-PRM can generate an extended demonstration based on the features extracted from the segmented original demonstration via 1D-SEG. This extended demonstration and the force item of the DMP model are fed into the BLS to update the DMP model to store previously learned skills (i.e., incremental learning). Three innovative characteristics of this study are represented below.

• Segmentation and extended demonstration: 1D-SEG is combined with G-PRM to generate an extended demonstration by incorporating the features of the original demonstration, so that fewer demonstrations are required to minimise the cost of data collection;
• DMP modelling and incremental learning update: The BLS learns the difference between the extended demonstration and the original demonstration by incrementally increasing the number of network nodes (hereafter referred to as additional enhancement nodes). The force item of the previously constructed DMP model is updated with the results generated by the BLS.
• Electric vehicle (EV) battery disassembly cases and an experimental platform were used to verify the effectiveness of the developed approach. Based on the approach, the successful disassembly of nuts and battery cells were achieved.

The remainder of this paper is structured as follows. Section 2 reviews the related work from the perspective of the trajectory generation and DMP-based optimisation. The key research gaps of the related work are summarised in Section 2. According to these gaps, segmentation and extended demonstration-based methods (namely 1D-SEG and G-PRM) and BLS-based incremental learning for DMP (namely Bi-DMP) are developed in Section 3. In Section 4, a disassembly experimentation platform based on EV batteries and collaborative robots is presented. The effectiveness of the proposed algorithms was verified by three cases, i.e., pick-and-place, unscrewing nuts, and battery cell removal. A comprehensive discussion is given in Section 5. Finally, Section 6 concludes the research.

2. Review of Related Work

2.1. Trajectory Generation

Trajectory generation is an important function to support robotic applications. Zhu et al. [3] developed an LfD-enabled trajectory generation algorithm that combines Gaussian noise with a task-parameterised Gaussian mixture model (TP-GMM). This algorithm facilitates the creation of the generation trajectory by augmenting the original demonstration. Hu et al. [21] proposed a method based on multi-region cost learning. By dividing a demonstration (trajectory) into various sub-regions and establishing cost functions for each sub-region, this method can fine-tune and preserve the geometric features of the demonstration more accurately to optimise the generated trajectory. Peng et al. [11] designed a generative adversarial imitation learning (GAIL)-based algorithm to optimise generated trajectories. In this algorithm, reinforcement learning and generative adversarial networks were combined to allow for more diversity in trajectory generation.

In order to achieve the local adjustment capability and meet the precision requirement of trajectory generation, Li et al. [12] proposed a trajectory generation method that combines a B-spline algorithm with a pruning algorithm. The pruning algorithm was used to create a trajectory, and B-spline was applied to smoothen the trajectory. However, the control points of the B-spline algorithm cannot accurately mimic the features of the created trajectory. To enhance the similarity of the generated trajectory to the original demonstration and enable local adjustment, sampling-based trajectory generation algorithms are highly promising choices, attributed to their simplicity, flexibility, and computational efficiency [13]. Hüppi et al. [13] and Zhou et al. [22] proposed a trajectory generation method based on PRM. PRM can generate trajectories through the Dijkstra algorithm and sampling points. The distribution of sampling points can affect both the algorithm’s efficiency and the quality of the generated trajectory. It is important to design a suitable sampling strategy.

However, several limitations persist when applying these approaches in manufacturing processes: (1) due to the lack of efficient trajectory segmentation methods, the generation of trajectories can only rely on the entire original demonstration, which may not be economical enough when trajectory adjustments are needed for just one sub-region. (2) Mimicking the characteristics of the demonstration proves challenging, and the shape of generated trajectories is hardly controlled because there is a lack of research on efficient sampling strategy and similarity measurement between the original demonstration and extended demonstration.

2.2. DMP-Based Optimisation

DMP has been widely used in LfD. Park et al. [23] and Ijspeert et al. [24] refined the classic DMP to accommodate demonstrations with the same start and target, as well as with obstacle avoidance, which provides a basis for further optimisation. The force item in DMP is constructed based on RBF kernels and weights which play a crucial role in determining the learning effect on demonstration. Consequently, the optimisation of the force item has received substantial attention. Teng et al. [25] integrated the Takagi–Sugeno fuzzy system with DMP, indicating that the fuzzy inference system can estimate the force item. To better learn the non-linear demonstrations, various researchers used neural networks to optimise the force item. For instance, Si et al. [15] adopted radial basis function neural networks (RBFNNs) to model the forcing item in DMP, enabling the learning of the position and orientation of demonstrations. To further improve the quality of DMP-generated trajectories, both Noohian et al. [16] and Kim et al. [17] not only replaced the force item with neural networks but also applied reinforcement learning for model training. The searching mechanisms in RL can improve the quality of trajectories generated by DMP. Moreover, Davchev et al. [6] introduced a residual correction policy aimed at improving the generalisation of DMP in peg-in-hole applications, thereby refining the quality of trajectories generated by DMP. Notably, the aforementioned methods are mainly focused on the improvement of force item and the DMP training process. Due to a lack of incremental learning capacity, when confronted with a new demonstration, the DMP has to be retrained.

To enable incremental learning within DMP, several scholars have made preliminary explorations. Lu et al. [18,19] attempted to combine a BLS with the force item of DMP, aiming to equip DMP with incremental learning capability. The BLS has a flat network structure and rapid update capability [20]. This allows the DMP to be iteratively updated efficiently. However, this approach is complex in force item construction. Simplifying the approach for combining BLS and DMP and reducing the computational complexity are necessary. In conclusion, optimising and enhancing the force item to endow DMP with incremental learning capabilities are crucial in real applications.

Therefore, the key research gaps of the related works can be summarised as follows:

• Mimicking the geometric features of the original demonstration proves challenging. There is a lack of research in designing an efficient sampling strategy and similarity measurement to avoid sharp turns in generated trajectories and enhance their resemblance to the demonstrations.
• These approaches neglect integrating the incremental learning function within DMP. When presented with a new demonstration, the DMP model has to be retrained. Meanwhile, there is a shortage of industrial applications leveraging incremental learning.

3. Research Methodology

The methodology in this research consists of the following steps: (i) segmentation and extended demonstration; (ii) DMP modelling and incremental learning update. The framework of this section is shown in Figure 2. Details are given below.

Before introducing the methodology of this paper, it is necessary to state the superiority of the designed method. Traditional trajectory segmentation methods mainly include subjective decision making based on trajectory geometric features [26,27] and adaptive segmentation based on neural networks [28]. In the optimised 1D-SEG developed in this paper, both the geometric features of the trajectory and a rigorous mathematical process are considered. This can make the segmentation results more objective. Compared to other trajectory generation algorithms like rapidly exploring random tree (RRT) [29], the advantage of PRM lies in the ease of optimising sampling points, as it constructs a roadmap to facilitate finding better trajectories in complex environments. In addition, the flat structure and rapid update strategy of the BLS do not require extensive training time, making it suitable for quickly building models.

3.1. Segmentation of Demonstration

Curvature serves as a quantitative measure of a trajectory’s features and, inspired by [30], a curvature-based trajectory segmentation algorithm has been designed. The details are as follows:

As illustrated in Figure 3, the blue dashed curve represents the curvatures of a trajectory. According to the changing trend of the curvatures, the trajectory can be divided into sub-regions represented as P_i. Ma_i and Mi_i denote the local maximum and local minimum curvatures in the ith region. Ma₀ and Mi₃ are the global maximum and global minimum curvature for the entire trajectory. The Euclidean distance ($Eudis$) between Ma_i and Mi_i in the ith region is defined as the curvature distance of the ith region ||P_i||. To optimise the segmentation of the trajectory, a threshold ε is defined. If ||P_i|| < ε, the ith region and its Ma_i and Mi_i will be eliminated in the segmentation process. The rest of the local maximum and local minimum curvatures are selected as the segmentation points.

3.2. Extended Demonstration via G-PRM

An extended demonstration will be created for the new target. The extended demonstration not only retains the features of the original demonstration but also reduces the computational complexity compared to traditional neural network algorithms. The process is detailed below.

After segmentation, the original demonstration $y_d$ is divided into several sub-demonstrations (specifically, sub-1 and sub-2 as illustrated in Figure 4b). Subsequently, the $Eudis$ from the segmentation point and the target point to the new target point are calculated as ||d₁|| and ||d₂||, respectively, as shown in Figure 4b. Since ||d₂|| < ||d₁||, sub-2 is closer to the new target point. Consequently, sub-2 is selected as the trajectory to be imitated (${Tra}_{imi}$) for generating the extended demonstration (${Tra}_{gen}$), as shown in Figure 4b,c.

For generating the ${Tra}_{gen}$, the PRM-based method is a preferred choice [31]. Figure 5 illustrates the fundamental processes of PRM. Initially, PRM involves generating sampling points (represented by blue dots) while avoiding obstacles, as depicted in Figure 5a. Subsequently, each sampling point is connected to its adjacent sampling points to form a roadmap, as illustrated in Figure 5b. Finally, to determine the shortest trajectory from position 1 to position 2, the Dijkstra algorithm is employed to search results in the roadmap. Figure 6 displays the flowchart of the G-PRM, with two red blocks highlighting the improvements. To more accurately express the meaning of the following equations, specific letters that represent sets and matrices will be presented in bold.

3.2.1. Generation of Sampling Points

According to Figure 5, the sampling points in PRM are generated randomly, which makes it difficult to control the shape of the resulting trajectory. By using Gaussian distribution-based sampling points (Gaussian sampling points, $\mathit{G}\mathit{S}\mathit{P}$), the adjustment of the generated trajectory ${Tra}_{gen}$ can be effectively facilitated.

The initialisation of $\mathit{G}\mathit{S}\mathit{P}$ relies on Gaussian noise and coordinate differences. Here are the details: a Gaussian distribution is defined as $N(\mu ,{\sigma }^2)$, where μ is the mean, and ${\sigma }^2$ represents the variance. Gaussian noise points directly generated from the $\mathcal{N}(\mu ,{\sigma }^2)$ are defined as $\mathit{G}\mathit{N}=\{{GN}_1,{GN}_2,\dots ,{GN}_{sample}\}$. The ‘$~$’ symbol indicates that each dimension of the element in the $\mathit{G}\mathit{N}$ follows $\mathcal{N}(\mu ,{\sigma }^2)$, which is expressed as Equation (1). $sample$ denotes the number of $\mathit{G}\mathit{N}$, and ${GN}_i=({GN}_{i_x},{GN}_{i_y},{GN}_{i_z})$, $i\in [1,sample]$.

$${GN}_{i_x}~\mathcal{N}\left(\mu ,{\sigma }^2\right),{GN}_{i_y}~\mathcal{N}\left(\mu ,{\sigma }^2\right),{GN}_{i_z}~\mathcal{N}\left(\mu ,{\sigma }^2\right)$$

(1)

$\mathit{G}\mathit{N}$ needs further processing to yield $\mathit{G}\mathit{S}\mathit{P}$, which is distributed between the segmentation point and the new target point. The number of elements in $\mathit{G}\mathit{S}\mathit{P}$ matches that in $\mathit{G}\mathit{N}$, hence $\mathit{G}\mathit{S}\mathit{P}=\{{GSP}_1, {GSP}_2,\dots ,{GSP}_{sample}\}$, and ${GSP}_i=({GSP}_{i_x},{GSP}_{i_y},{GSP}_{i_z})$, with the same $i\in [1,sample]$. A time series $\mathit{t}=\{t_1,t_2,\dots ,t_{sample}\}$ serves as a control variable to calculate the ${GSP}_i$ from ${GN}_i$ based on the coordinate differences between the segmentation point and the new target point. The segmentation point and new target point for ${Tra}_{gen}$ are $seg=\left(x_s{,y}_s,z_s\right)$ and $n\_tar=\left(x_e{,y}_e,z_e\right)$, respectively. The calculation for $t_i$ in $\mathit{t}$ and ${GSP}_i$ is described in Equation (2):

$$\left\{\begin{array}{c}{GSP}_{i_x}=x_s+t_i·\left(x_e-x_s\right)+{GN}_{i_x}\\ {GSP}_{i_y}=y_s+t_i·\left(y_e-y_s\right)+{GN}_{i_y}\\ G{SP}_{i_z}=z_s+t_i·\left(z_e-z_s\right)+{GN}_{i_z}\end{array}\right., t_i=\frac{i}{sample}, i=0,1,\dots ,sample$$

(2)

3.2.2. Bias Optimisation

G-PRM can generate an initial ${Tra}_{gen}$ based on $\mathit{G}\mathit{S}\mathit{P}$. To derive bias $\mathit{b}\mathit{v}$ for updating $\mathit{G}\mathit{S}\mathit{P}$ and ${Tra}_{gen}$, two key processes are required: ➀ it is essential to design a Fréchet distance-based [32] similarity criterion to measure the similarity between ${Tra}_{gen}$ and ${Tra}_{imi}$ ($Fresim({Tra}_{gen},{Tra}_{imi})$). $Fresim({Tra}_{gen},{Tra}_{imi})$ is derived from Fréchet distance. ➁ Scaling ${Tra}_{gen}$ and ${Tra}_{imi}$ to the same scale facilitates the calculation of coordinate differences.

➀ To calculate the Fréchet distance between ${Tra}_{imi}$ and ${Tra}_{gen}$, three steps are involved:

Step 1—Initialisation definitions: This paper defines $Fré\left(a,b\right)$ as the Fréchet distance between the first $a$ points of ${Tra}_{imi}$ and the first $b$ points of ${Tra}_{gen}$. Both ${Tra}_{imi}$ and ${Tra}_{gen}$ contain an equal number of points, with $a$ and $b\in (1, k]$. $k$ represents the number of points on the two trajectories, respectively. Since ${Tra}_{imi}$ and ${Tra}_{gen}$ are discrete, the computation of the Fréchet distance of the two trajectories is based on $Eudis$.

Step 2—Boundary conditions for Fréchet distance:

$$\begin{array}{l}Fr\mathrm{é}\left(1,1\right)=Eudis\left(seg,seg\right)\\ Fr\mathrm{é}\left(a,1\right)=max\left(Fr\mathrm{é}\left(a-1,1\right),Eudis\left({Tra}_{imi\_a},seg\right)\right),a>1\\ Fr\mathrm{é}\left(1,b\right)=max\left(Fr\mathrm{é}\left(1,b-1\right),Eudis\left(seg,{Tra}_{gen\_b}\right)\right),b>1\end{array}$$

(3)

where $Fré\left(1,1\right)$ in Equation (3) represents the Fréchet distance between the first points ($seg$) of ${Tra}_{imi}$ and ${Tra}_{gen}$, the $Eudis\left(·\right)$ represents the distance between two elements in the bracket. ${Tra}_{imi\_a}$ and ${Tra}_{gen\_b}$ represent the ath point and the bth point in ${Tra}_{imi}$ and ${Tra}_{gen}$, respectively. The $\mathit{max}(·)$ the largest element in the brackets.

Step 3—Recursive calculation for Fréchet distance: Based on Equation (3), the Fréchet distance between ${Tra}_{imi}$ and ${Tra}_{gen}$ can be calculated as shown in Equation (4):

$$Fr\mathrm{é}\left(a,b\right)=max\left(\begin{array}{c}\min \left(Fr\mathrm{é}\left(a-1,b\right),Fr\mathrm{é}\left(a,b-1\right),Fr\mathrm{é}\left(a-1,b-1\right)\right)\\ , Eudis\left({Tra}_{imi\_a},{Tra}_{gen\_b}\right)\end{array}\right),a>1 and b>1$$

(4)

where $\min (·)$ opts for the smallest value in brackets. When $a=k$ and $b=k$, the Fréchet distance $Fré\left({Tra}_{imi},{Tra}_{gen}\right)$ between ${Tra}_{imi}$ and ${Tra}_{gen}$ is obtained.

Based on the calculation results from Equation (4), $Fresim({Tra}_{gen},{Tra}_{imi})$ is calculated as shown in Equation (5):

$$Fresim({Tra}_{gen},{Tra}_{imi})=max(1-\frac{Fr\mathrm{é}({Tra}_{gen},{Tra}_{imi})}{\frac{\sqrt{len({Tra}_{gen}) \times len({Tra}_{imi})}}{\sqrt{2}}},0)$$

(5)

where $\sqrt{len\left({Tra}_{gen}\right)\times len\left({Tra}_{imi}\right)}$ is the geometric mean of the lengths of the ${Tra}_{imi}$ and ${Tra}_{gen}$. $len(·)$ represents the length of the trajectory calculated via $Eudis$. $Fresim(·)$ ∈ [0, 1]. A score of 0 indicates complete dissimilarity between the two trajectories, whereas 1 signifies that the trajectories are identical.

➁ To facilitate bias optimisation, ${Tra}_{gen}$ needs to be scaled to the same scale as ${Tra}_{imi}$. The target point of ${Tra}_{imi}$ is represented by $tar=\left(x_o{,y}_o,z_o\right)$ as shown by the red dot in Figure 4. The calculated scaling vector is denoted as $\overrightarrow{scal{e}_{factor}}$, as shown in Equation (6):

$$\overrightarrow{scal{e}_{factor}}=\overrightarrow{(\frac{{x_e-x}_s}{{x_o-x}_s},\frac{{y_e-y}_s}{{y_o-y}_s},\frac{{z_e-z}_s}{{z_o-z}_s})}$$

(6)

$\overrightarrow{scal{e}_{factor}}$ includes scaling factors for x, y, and z axes, denoted as $scal{e}_{factor\_x}$, $scal{e}_{factor\_y}$, and $scal{e}_{factor\_z}$, respectively. The scaled trajectory of ${Tra}_{gen}$ is defined as ${Tra}_s$. Each coordinate in ${Tra}_s$ can be calculated by Equation (7):

$$\left\{\begin{array}{c}{Tra}_{{s\_i}_x}=x_s+scal{e}_{factor\_x}·({Tra}_{{gen\_i}_x}-x_s)\\ {Tra}_{{s\_i}_y}=y_s+scal{e}_{factor\_y}·({Tra}_{{gen\_i}_y}-y_s)\\ {Tra}_{{s\_i}_z}=z_s+scal{e}_{factor\_z}·({Tra}_{{gen\_i}_z}-z_s)\end{array}\right.,i=(1,2\dots ,k)$$

(7)

where ${Tra}_{{s\_i}_x}$, ${Tra}_{{s\_i}_y}$, and ${Tra}_{{s\_i}_z}$ represent the x, y, and z coordinates of the ith point in ${Tra}_s$. ${Tra}_{{gen\_i}_x}$, ${Tra}_{{gen\_i}_y},$ and ${Tra}_{{gen\_i}_z}$ refer to the x, y, and z coordinates of the ith point in ${Tra}_{gen}$.

Finally, the bias $\mathit{b}\mathit{v}=\{{bv}_1,{bv}_2,\dots ,{bv}_k\}$ for updating $\mathit{G}\mathit{S}\mathit{P}$ is calculated based on the results of ➀ and ➁. The initial $\mathit{b}\mathit{v}$ in each dimension is zero. $\rho$ is the learning factor, and $iter$ represents the current number of iterations. Equation (8) shows the calculation process of $\mathit{b}\mathit{v}$:

$$\left\{\begin{array}{c}{bv}_{i_x}\leftarrow {bv}_{i_x}+{(1-Fresim)·\rho }^{iter}·({Tra}_{{imi\_i}_x}-{Tra}_{{s\_i}_x})\\ {bv}_{i_y}\leftarrow {bv}_{i_y}+{(1-Fresim)·\rho }^{iter}·\left({Tra}_{{imi\_i}_y}{-Tra}_{{s\_i}_y}\right) \\ {bv}_{i_z}\leftarrow {bv}_{i_z}+{(1-Fresim)·\rho }^{iter}·\left({Tra}_{{imi\_i}_z}{-Tra}_{{s\_i}_z}\right)\end{array}\right.,i=(1,2\dots ,k)$$

(8)

In Equation (8), ${bv}_{i_x}$, ${bv}_{i_y}$, and ${bv}_{i_z}$ represent coordinate differences between ${Tra}_s$ and ${Tra}_{imi}$ in the x, y, and z dimensions, respectively. ${Tra}_{{imi\_i}_x}$, ${Tra}_{{imi\_i}_y}$, and ${Tra}_{{imi\_i}_z}$ are the x, y, and z coordinates of the ith point in ${Tra}_{imi}$. To update $\mathit{G}\mathit{S}\mathit{P}$ based on the $\mathit{b}\mathit{v}$, interpolation is a necessary method to align the number of elements in $\mathit{b}\mathit{v}$ with $\mathit{G}\mathit{S}\mathit{P}$. After interpolation, the $\mathit{b}\mathit{v}$ is renamed and reformulated as $\mathit{I}\mathit{b}\mathit{v}=\{{Ibv}_1,{Ibv}_2,\dots ,{Ibv}_{sample}\}$. Subsequently, updated $\mathit{G}\mathit{S}\mathit{P}$ ($\mathit{G}\mathit{S}\mathit{P}\mathit{n}\mathit{e}\mathit{w}$) is generated through Equation (9), where each coordinate is calculated accordingly:

$$\left\{\begin{array}{c}{GSPnew}_{i_x}={GSP}_{i_x}+{Ibv}_{i_x}\\ G{SPnew}_{i_y}={GSP}_{i_y}+{Ibv}_{i_y}\\ {GSPnew}_{i_z}={GSP}_{i_z}+{Ibv}_{i_z}\end{array}\right.,i=0,1,\dots ,sample$$

(9)

The $\mathit{G}\mathit{S}\mathit{P}\mathit{n}\mathit{e}\mathit{w}$ is used to generate updated ${Tra}_{gen}$ via G-PRM. $\mathit{G}\mathit{S}\mathit{P}$ is regenerated based on Equations (1) and (2) during each iteration. Subsequently, $Fresim(·)$ needs to be recalculated based on updated ${Tra}_{gen}$ and ${Tra}_{imi}$. Additionally, Equations (6)–(9) are applied iteratively to continue updating ${Tra}_{gen}$ and $\mathit{G}\mathit{S}\mathit{P}\mathit{n}\mathit{e}\mathit{w}$. Through this process, the similarity between ${Tra}_{gen}$ and ${Tra}_{imi}$ gradually increases. After reaching the pre-set number of iterations, bias optimisation is completed.

3.3. DMP Modelling and Incremental Learning Updating

3.3.1. Modelling of DMP

In this session, the $y_d$ is used as an example to model the DMP. This study focuses on discrete DMP represented in the Cartesian space. The Cartesian space representation exhibits advantages over the joint space representation as it allows for trajectory planning without considering the joints and their relative positions [33]. The DMP model is represented below [4,23,24]:

$${\mathrm{\tau }}^2{\ddot{y}}_d={\alpha }_{y_d}\left({\beta }_{y_d}(g-y_d)-\tau {\dot{y}}_d\right)+Kf(\vartheta )-K\vartheta (g-y_{d0})$$

(10)

where $y_d$ represents the original demonstration; ${\dot{y}}_d$ and ${\ddot{y}}_d$ denote the first- and second-order derivatives of the demonstration, respectively; ${\alpha }_{y_d}$ and ${\beta }_{y_d}$ (${{\alpha }_{y_d}, \beta }_{y_d}>0$) are damping and stiffness factors (normally ${\alpha }_{y_d}=4{\beta }_{y_d}$); $\tau$ is the timing parameter used to adjust the timestep of the demonstration; $y_{d0}$ and $g$ mark the start and target of the demonstration; $K$ is a constant that is equal to ${\alpha }_{y_d}$; $\vartheta$ is the phase variable computed by the canonical system as expressed using Equation (11). It is used to calculate the value of the force item for each timestep; $f\left(\vartheta \right)$ is the force item represented in Equation (12); ${\omega }_i$ in $f\left(\vartheta \right)$ is calculated by locally weighted regression (LWR) [24].

$$\tau \dot{\vartheta }=-{\alpha }_{\vartheta }\vartheta$$

(11)

$$f(\vartheta )=\frac{{\sum }_{i=1}^N{\omega }_i{\mathsf{\psi }}_i\left(\vartheta \right)}{{\sum }_{i=1}^N{\mathsf{\psi }}_i(\vartheta )}\vartheta , {\mathsf{\psi }}_i(\vartheta )=\exp (-h_i\left(\vartheta -c_i{)}^2\right)$$

(12)

In Equations (11) and (12), ${\alpha }_{\vartheta }$ is a constant, and it is always set to one; $\tau$ is the same as that defined in Equation (10); $N$ denotes the number of RBF kernels; ${\mathsf{\psi }}_i\left(\vartheta \right)$ represents the ith RBF kernel; $c_i>0$ and $h_i>0$ are the centre and width of the ith RBF kernel.

For computational convenience, Equation (11) is transformed into Equation (13) to calculate $\vartheta$:

$$\begin{array}{c}d\vartheta =-\tau {\alpha }_{\vartheta }\vartheta dt\\ \vartheta =\vartheta +d\vartheta \end{array}$$

(13)

where $dt$ is the inverse of the number of the demonstration’s timestep points. When designating a new target $g_{new}$ to replace $g$ in Equation (10), each ${\ddot{y}}_d$ in the generated trajectory, denoted as ${\ddot{y}}_{{gen}_l}$, is calculated by using Equation (10). It is updated with the change in $\vartheta$; $l$ is a timestep within the generated trajectory.

The next coordinate $y_{{gen}_{l+1}}$ is calculated by using Equation (14):

$${\dot{y}}_{{gen}_{l+1}}\leftarrow {\dot{y}}_{{gen}_l}+{\tau \ddot{y}}_{{gen}_l}dt, y_{{gen}_{l+1}}\leftarrow y_{{gen}_l}+{\tau \dot{y}}_{{gen}_l}dt$$

(14)

3.3.2. Bi-DMP for Incremental Learning Updating

In this section, the structure and mathematical formulation of Bi-DMP are constructed. The core task of Bi-DMP is to construct the incremental force item for DMP. Inspired by Lu et al. [18,19] and Chen et al. [20], the main structure of the BLS designed for DMP is shown in Figure 7.

In Figure 7, the grey dots symbolise the output of the BLS. The hidden layer of the BLS is represented by blue, red, and orange dots. The blue dots denote mapping features, while the red and orange dots correspond to enhancement nodes and additional enhancement nodes, respectively. Therefore, the incremental force item is the hidden layer combined with the weight corresponding to $\mathbf{∆}{\mathit{F}}_{\mathit{I}\mathit{n}}$. Further explanations of the remaining labels in Figure 7 will be provided in the following sections.

This paper postulates that the force item of the $y_d$ is ${\mathbf{F}}_{\mathit{o}\mathit{r}\mathit{i}}$. ${\mathbf{F}}_{\mathit{o}\mathit{r}\mathit{i}}$ has been computed and parameterised by LWR [24] and is represented in Equation (12). The ${Tra}_{gen}$ should combined with sub-1 to form a complete trajectory (namely ${Tra}_{gen}$-based complete demonstration) as shown in Figure 4c, namely $y_{dn}$. It is used to calculate the original target force item ${\mathbf{F}}_{\mathit{t}\mathit{a}\mathit{r}\mathit{g}\mathit{e}\mathit{t}}$. Both ${\mathbf{F}}_{\mathit{o}\mathit{r}\mathit{i}}$ and ${\mathbf{F}}_{\mathit{t}\mathit{a}\mathit{r}\mathit{g}\mathit{e}\mathit{t}}$ are sets containing the force item value for each point in demonstrations. Replacing $y_d$ with $y_{dn}$ in Equation (10), the initial ${\mathbf{∆}\mathbf{F}}_{\mathit{I}\mathit{n}}$ can be written as Equation (15):

$$\begin{array}{l}\mathbf{∆}{\mathit{F}}_{\mathit{I}\mathit{n}}={\mathit{F}}_{\mathit{t}\mathit{a}\mathit{r}\mathit{g}\mathit{e}\mathit{t}}-{\mathit{F}}_{\mathit{o}\mathit{r}\mathit{i}}\\ =\frac{{\mathrm{\tau }}^2{\ddot{y}}_{dn}-{\alpha }_{y_{dn}}\left({\beta }_{y_{dn}}\left(g_{y_{dn}}-y_{dn}\right)-\tau {\dot{y}}_{dn}\right)}{K}+\vartheta \left(g_{y_{dn}}-y_{dn0}\right)-f(\vartheta )\end{array}$$

(15)

The key components of the hidden layer are ${\mathit{Z}}_{\mathit{i}}$, $\mathit{H}$, and $\mathit{E}$, respectively. ${\mathit{Z}}_{\mathit{i}}$ refers to the mapping features, $\mathit{H}$ represents the enhancement nodes, and $\mathit{E}$ indicates the additional enhancement nodes. This paper assumes that $u$ is the number of mapping features, $l$ denotes the number of enhancement nodes, and $s$ indicates the number of additional enhancement nodes, as depicted from Equation (16) to Equation (18):

$${\mathit{Z}}_{\mathit{i}}=\delta \left(\mathit{\phi }{\left(\vartheta \right)\mathit{W}}_{{\mathit{Z}}_{\mathit{i}}}+{\mathit{\beta }}_{{\mathit{Z}}_{\mathit{i}}}\right),i=1,2,\dots ,u$$

(16)

$$\mathit{H}=\xi \left(\mathbf{\Xi }(\vartheta ){\mathit{W}}_{{\mathit{H}}_{\mathit{l}}}+{\mathit{\beta }}_{{\mathit{H}}_{\mathit{l}}}\right)$$

(17)

$$\mathit{E}=\xi \left(\mathbf{\Xi }(\vartheta ){\mathit{W}}_{{\mathit{E}}_{\mathit{s}}}+{\mathit{\beta }}_{{\mathit{E}}_{\mathit{s}}}\right)$$

(18)

${\mathit{W}}_{{\mathit{Z}}_{\mathit{i}}}$ is a randomly initialised weight matrix and subsequently fine-tuned via a sparse autoencoder [20]. Both ${\mathit{W}}_{{\mathit{H}}_{\mathit{l}}}$ and ${\mathit{W}}_{{\mathit{E}}_{\mathit{s}}}$ are initialised as orthogonal matrixes. The terms ${\mathit{\beta }}_{{\mathit{Z}}_{\mathit{i}}}$, ${\mathit{\beta }}_{{\mathit{H}}_{\mathit{l}}}$, and ${\mathit{\beta }}_{{\mathit{E}}_{\mathit{s}}}$ represent biases, while $\delta$ and $\xi$ are activation functions. As suggested in [20], $\delta$ can be a linear activation function, while $\xi$ can be the Tansig function, as detailed in Equation (19):

$$\xi \left(·\right)=\frac{2}{1+e^{-2\left(·\right)}}-1$$

(19)

To determine the weight ${\mathit{W}}^{\mathbf{∆}{\mathit{F}}_{\mathit{I}\mathit{n}}}$ corresponding to $\mathbf{∆}{\mathit{F}}_{\mathit{I}\mathit{n}}$, it is necessary to calculate the basis weight matrix ${\mathit{W}}^{\mathit{b}}$ based on mapping features and enhancement nodes. Equations (20)–(22) represent the computation process:

$$\begin{array}{l}\mathbf{∆}{\mathit{F}}_{\mathit{I}\mathit{n}}=\left[{\mathit{Z}}_{\mathbf{1}},\dots ,{\mathit{Z}}_{\mathit{u}}|\xi \left(\mathbf{\Xi }\left(\vartheta \right){\mathit{W}}_{{\mathit{H}}_{\mathit{l}}}+{\mathit{\beta }}_{{\mathit{H}}_{\mathit{l}}}\right)\right]{\mathit{W}}^{\mathit{b}}\\ =\left[\mathit{Z}|\mathit{H}\right]{\mathit{W}}^{\mathit{b}}\end{array}$$

(20)

$${\mathit{W}}^{\mathit{b}}={\left[\mathit{Z}|\mathit{H}\right]}^{-1}\mathbf{∆}{\mathit{F}}_{\mathit{I}\mathit{n}}={\mathcal{A}}^{\mathfrak{f}}\mathbf{∆}{\mathit{F}}_{\mathit{I}\mathit{n}}$$

(21)

$${\mathcal{A}}^{\mathfrak{f}}=\underset{\lambda \to 0}{\lim }{(\lambda \mathit{I}+\mathcal{A}{\mathcal{A}}^{\mathit{T}})}^{-1}{\mathcal{A}}^{\mathit{T}}$$

(22)

In Equation (21), $\mathcal{A}$ denotes matrix $\left[\mathit{Z}|\mathit{H}\right]$, where ${\mathcal{A}}^{\mathfrak{f}}$ represents the Moore–Penrose pseudo-inverse matrix of $\mathcal{A}$. ${\mathcal{A}}^{\mathit{T}}$ is a transpose matrix of $\mathcal{A}$, and $\lambda$ is the regularisation parameter. I is an identity matrix. Once the ${\mathit{W}}^{\mathit{b}}$ is obtained, the $\mathit{E}$ is incorporated into Equation (20) to improve ${\mathit{W}}^{\mathit{b}}$ for higher accuracy in approximating $\mathbf{∆}{\mathit{F}}_{\mathit{I}\mathit{n}}$. Then, Equation (20) is rewritten as Equation (23):

$$\begin{array}{l}∆{\mathit{F}}_{\mathit{I}\mathit{n}}=\left[\mathit{Z}|\mathit{H}|\xi \left(\mathbf{\Xi }\left(\vartheta \right){\mathit{W}}_{{\mathit{E}}_{\mathit{s}}}+{\mathit{\beta }}_{{\mathit{E}}_{\mathit{s}}}\right)\right]{\mathit{W}}^{∆{\mathit{F}}_{\mathit{I}\mathit{n}}}\\ \begin{array}{c} =\left[\mathcal{A}|\mathit{E}\right]{\mathit{W}}^{∆{\mathit{F}}_{\mathit{I}\mathit{n}}}\\ =\mathcal{A}\mathit{n}\mathit{e}\mathit{w}{\mathit{W}}^{∆{\mathit{F}}_{\mathit{I}\mathit{n}}}\end{array}\end{array}$$

(23)

where $\mathcal{A}\mathit{n}\mathit{e}\mathit{w}=\left[\mathcal{A}|\mathit{E}\right]$ (hidden layer matrix), The advantage of incremental learning is that when $\mathit{E}$ is added, there is no need to recalculate ${\left[\mathcal{A}\right|\mathit{E}]}^{-1}$, thereby accelerating the solution process and facilitating network expansion [20]. The ${\left[\mathcal{A}\right|\mathit{E}]}^{-1}$ is defined as ${\mathcal{A}\mathit{n}\mathit{e}\mathit{w}}^{\mathfrak{f}}$. ${\mathcal{A}\mathit{n}\mathit{e}\mathit{w}}^{\mathfrak{f}}$ is calculated by Equation (24):

$${\mathcal{A}\mathit{n}\mathit{e}\mathit{w}}^{\mathfrak{f}}=[\begin{array}{c}{\mathcal{A}}^{\mathfrak{f}}-\mathit{D}{\mathit{B}}^{\mathit{T}}\\ {\mathit{B}}^{\mathit{T}}\end{array}]$$

(24)

The temporary matrix $\mathit{B}$ and $\mathit{D}$ can be calculated according to Equation (25):

$$\mathit{B}=\left\{\begin{array}{c}{\mathit{C}}^{\mathfrak{f}} , if \mathit{C}\ne 0\hfill \\ {(\mathit{I}+{\mathit{D}}^{\mathit{T}}\mathit{D})}^{-1}{\mathit{D}}^{\mathit{T}}{\mathcal{A}}^{\mathfrak{f}}, if \mathit{C}=0\end{array}\right., \mathit{D}={\mathcal{A}}^{\mathfrak{f}}\mathit{E}$$

(25)

where $\mathit{C}=\mathit{E}-\mathcal{A}\mathit{D}$. ${\mathit{C}}^{\mathfrak{f}}$ represents the Moore–Penrose pseudo-inverse matrix of $\mathit{C}$. Therefore, the ${\mathit{W}}^{\mathbf{∆}{\mathit{F}}_{\mathit{I}\mathit{n}}}$ can be calculated according to Equation (23) as Equation (25):

$${\mathit{W}}^{∆{\mathit{F}}_{\mathit{I}\mathit{n}}}{=\mathcal{A}\mathit{n}\mathit{e}\mathit{w}}^{\mathfrak{f}}∆{\mathit{F}}_{\mathit{I}\mathit{n}}=[\begin{array}{c}{\mathit{W}}^{\mathit{b}}-\mathit{D}{\mathit{B}}^{\mathit{T}}∆{\mathit{F}}_{\mathit{I}\mathit{n}}\\ {\mathit{B}}^{\mathit{T}}∆{\mathit{F}}_{\mathit{I}\mathit{n}}\end{array}]$$

(26)

3.3.3. Iterative Update of Force Item under Multi-Demonstrations

Finally, the force item $f\left(\vartheta \right)$ can be rewritten as incremental force item $f_{in}\left(\vartheta \right)$. According to Equations (16)–(18), $\vartheta$ in DMP can modify each element within the BLS. Consequently, to maintain consistency in $f_{in}\left(\vartheta \right)$, $\mathcal{A}\mathit{n}\mathit{e}\mathit{w}$ is reformulated as $\mathcal{A}\mathit{n}\mathit{e}\mathit{w}\left(\vartheta \right)=\left[{\mathit{Z}}_{\mathbf{1}},\dots ,{\mathit{Z}}_{\mathit{u}}|\xi \left(\mathbf{\Xi }\left(\vartheta \right){\mathit{W}}_{{\mathit{H}}_{\mathit{l}}}+{\mathit{\beta }}_{{\mathit{H}}_{\mathit{l}}}\right)|\xi \left(\mathbf{\Xi }\left(\vartheta \right){\mathit{W}}_{{\mathit{E}}_{\mathit{s}}}+{\mathit{\beta }}_{{\mathit{E}}_{\mathit{s}}}\right)\right]$. Equation (27) shows the $f_{in}\left(\vartheta \right)$.

$$\begin{array}{l}{f}_{in}(\vartheta )=\frac{{\sum }_{i=1}^N{\vartheta \mathsf{\psi }}_i\left(\vartheta \right)}{{\sum }_{i=1}^N{\mathsf{\psi }}_i(\vartheta )}{\omega }_i+\mathcal{A}\mathit{n}\mathit{e}\mathit{w}(\vartheta ){\mathit{W}}^{\mathbf{∆}{\mathit{F}}_{\mathit{I}\mathit{n}}}\\ =\mathit{\phi }\left(\vartheta \right)\mathit{W}+\mathcal{A}\mathit{n}\mathit{e}\mathit{w}(\vartheta ){\mathit{W}}^{\mathbf{∆}{\mathit{F}}_{\mathit{I}\mathit{n}}}\end{array}$$

(27)

When there are additional demonstrations, $f_{in}\left(\vartheta \right)$ should be revised to ${f_{in}}_{\mathcal{P}}\left(\vartheta \right)$, where $\mathcal{P}$ is the number of demonstrations. Additionally, the $\mathbf{∆}{\mathit{F}}_{\mathit{I}\mathit{n}}$ associated with the new demonstration equals the difference between the current $f_{in}\left(\vartheta \right)$ and the new demonstration’s ${\mathit{F}}_{\mathit{t}\mathit{a}\mathit{r}\mathit{g}\mathit{e}\mathit{t}}$. The new difference can be calculated based on Equation (15). Then, the new ${\mathit{W}}^{\mathbf{∆}{\mathit{F}}_{\mathit{I}\mathit{n}}}$ (namely ${\mathit{W}}_{\mathcal{p}}^{\mathbf{∆}{\mathit{F}}_{\mathit{I}\mathit{n}}}, \mathcal{p}=1,2,\dots \mathcal{P}$) and new $\mathcal{A}\mathit{n}\mathit{e}\mathit{w}\left(\vartheta \right)$ (namely ${\mathcal{A}\mathit{n}\mathit{e}\mathit{w}\left(\vartheta \right)}_{\mathcal{p}}, \mathcal{p}=1,2,\dots \mathcal{P}$) are determined according to Equations (20)–(26) based on the new $\mathbf{∆}{\mathit{F}}_{\mathit{I}\mathit{n}}$. Consequently, the ${f_{in}}_{\mathcal{P}}\left(\vartheta \right)$ can be expressed as Equation (28). Moreover, when $\mathcal{P}=1$, Equation (28) simplifies to Equation (27).

$${f_{in}}_{\mathcal{P}}(\vartheta )=\mathit{\phi }\left(\vartheta \right)\mathit{W}+{\sum }_{\mathcal{p}=1}^{\mathcal{P}}{\mathcal{A}\mathit{n}\mathit{e}\mathit{w}\left(\vartheta \right)}_{\mathcal{p}}{\mathit{W}}_{\mathcal{p}}^{\mathbf{∆}{\mathit{F}}_{\mathit{I}\mathit{n}}},(\mathcal{p}=1,2,\dots \mathcal{P})$$

(28)

4. Experiments and Case Studies

In simulations, trajectories from the LASA dataset [34] and authors’ handwriting trajectories were adopted to validate the algorithm.

4.1. Extended Demonstration Based on 1D-SEG and G-PRM

Handwriting A-shaped trajectories and M/S-shaped trajectories from LASA are selected to validate the 1D-SEG and G-PRM. Initially, 1D-SEG is applied to segment the original demonstrations based on their curvatures, as illustrated from Figure 8a,b and Figure 9a. Subsequently, ${Tra}_{imi}$ is calculated according to 3.1, depicted by the purple dashed curves in Figure 9b. Following the algorithm in 3.2, G-PRM generates ${Tra}_{gen}$ (yellow curves in Figure 9b). More specifically, the red and green dots in Figure 9b represent the segmentation and new target points of ${Tra}_{gen}$, respectively. The blue dots are the target points of the original demonstrations.

4.2. Modelling of DMP

In this subsection, three trajectories are selected to construct DMP. As illustrated in Figure 10, the blue curve trajectories represent demonstrations, whereas the red dashed curve trajectories are the outcomes of DMP. The force item of the DMP will serve as a benchmark for the following incremental learning processes.

4.3. Bi-DMP for Extended Demonstration

Figure 11 intuitively expresses the geometric differences between original and extended demonstrations and the learning performance of Bi-DMP. The red dashed curves in Figure 11 represent the trajectories generated via Bi-DMP.

4.4. Case Study 1—Pick-and-Place-Based Bi-DMP

As the most common industrial step, pick-and-place is always used as a benchmark trajectory to validate algorithms. In this part of the study, operators of different ages and genders were selected to perform pick-and-place experiments to collect their demonstration trajectories. Specifically, the operators were asked to use a gripper to catch an object, drag the robot to move, and place the object into an orange box. Figure 12 shows the details:

Figure 12a shows the experimental platform, while Figure 12b presents the demonstrations from operators. These demonstrations were collected via the TM collaborative robot, and the trajectories are shown in Figure 13a. It can be observed that there are discrepancies in the demonstration trajectories when performing the same task, due to factors such as gender, height, and arm span [35,36]. These factors should not be ignored, as operators prefer to work with a robot that matches their habits and physical characteristics in HRC. For example, the trajectory demonstrated by a person of short stature may collide with the body of a person of tall stature. Therefore, using only one DMP to generalise this pick-and-place task is not appropriate. Bi-DMP, on the other hand, can model and retain the differences in all demonstrations.

In the experiment, the red curve trajectory is chosen as a benchmark demonstration for the construction of DMP, as shown in Figure 13b. The purple dashed curve in Figure 13c is the path generated by the DMP. In order to preserve the discrepancies between different operators, Bi-DMP is used to learn the two new trajectories based on the constructed DMP as shown in Figure 13d. The orange and blue dashed curves are generated by Bi-DMP.

According to the above results, Bi-DMP can effectively extend the force item of DMP. Starting from the constructed DMP, the force item is continuously updated so that the differences in the demonstrations can be fully learned.

4.5. Case Study 2—Unscrewing Nuts and Battery Cell Disassembly

4.5.1. Experiment Platform and Problem Descriptions

Unscrewing nuts and removing battery cells are critical steps in EV battery disassembly applications. The experimental platform is shown in Figure 14a. A multi-purpose adaptor and a vacuum chuck are designed and employed to perform the mentioned tasks, as shown in Figure 14b and Figure 14c, respectively.

Due to the hazardous nature of the EV battery and the constraints of the working space, it is imperative for operators to minimise the duration of the demonstration. This results in nuts closer to the original demonstration being easily disassembled, as illustrated by the red curves in Figure 14d. However, nuts positioned far from the demonstration bring challenges for disassembly due to generalisation issues, as illustrated by the red dashed trajectory in Figure 14d. In the case of battery cell removal, trajectories generated from the Bi-DMP previously used for unscrewing nuts exhibit a deviation due to the lack of specific demonstrations which may cause unforeseen collisions with the working environment. The algorithms proposed in this paper can address these problems.

4.5.2. Design and Use of End-Effectors

Figure 15a,b show the structural design and application of end-effectors for unscrewing nuts and removing battery cells, respectively. The end-effector designed for unscrewing nuts is a multi-purpose adaptor, comprising a series of steel cylinders and a spring system, as depicted in Figure 15a. During the unscrewing process, the steel cylinders adaptively enclose the nuts, allowing for their removal through rotation and friction induced by the multi-purpose adaptor. The end-effector for battery cell removal is a vacuum chuck. It consists of four chucks and pipelines, as shown in Figure 15b. To remove a battery cell, two chucks secure it, as highlighted in the red rectangular box in Figure 15b. Due to safety and payload considerations, only one battery cell should be removed at a time.

4.5.3. Demonstration and Implementation

To ensure the safety of the operator, only the nut at position No. 3 was unscrewed during the demonstration, as shown in Figure 16a. The demonstration adhered strictly to the disassembly rules: moving the multi-purpose adaptor from the start point to the target position (nut), unscrewing it, and then returning to the start place.

The diameter of the multi-purpose adaptor is 25 mm, and the length of the groove is 28.5 mm, as shown in Figure 16c. If the distance between the central area of the nut and the endpoint of the DMP-generated trajectory exceeds the threshold value of 1.75 mm, the unscrewing nuts task is likely to fail due to collisional interference, as highlighted in the red ellipse in Figure 16c. Figure 16b illustrates the generated trajectory via DMP. Positions No. 1 to No. 8 corresponding to the light blue trajectories generated by DMP were able to complete the unscrewing nuts task. In contrast, positions from No. 9 to No. 14 (purple trajectories) failed to complete the task due to excessive errors. Thus, the EV battery is categorised as having a successful area and a failed area.

Table 1. 1. Error between the central area of nut and the target point generated by DMP. 2. Error between the central area of nut and the target point generated by DMP (Continued).

1
Error (mm)	No. 1	No. 2	No. 3	No. 4	No. 5	No. 6	No. 7
	0.404	0.260	0	0.989	1.023	1.069	1.487
2
Error (mm)	No. 8	No. 9	No. 10	No. 11	No. 12	No. 13	No. 14
	1.431	1.751	1.760	1.787	2.760	2.743	2.730

shows the error associated with each position.

In this scenario, 1D-SEG and G-RPM are used to segment the demonstration and generate ${Tra}_{gen}$. The outcomes of 1D-SEG are shown in Figure 17a. Typically, the new target point for ${Tra}_{gen}$ is selected to be at the centre of the failed nut positions, marked by the brown pentagram in Figure 17b. This position provides better coverage of the working area. The yellow trajectory in Figure 17b shows the ${Tra}_{gen}$-based complete demonstration. Subsequently, Bi-DMP is trained based on the force item of previously constructed DMP. The trajectories generated by Bi-DMP are used to re-execute the unscrewing nuts task in the failed area as shown in Figure 17c. According to the error analysis provided in

Table 2. Error between the central area of nut and the target point generated by Bi-DMP.

Error (mm)	No. 9	No. 10	No. 11	No. 12	No. 13	No. 14
	1.477	1.248	1.477	1.477	1.248	1.477

, all errors are smaller than the threshold. Therefore, the remaining tasks can be executed successfully. Figure 18 shows the implementation of unscrewing nuts.

Another significant application is the removal of battery cells. However, due to the height difference between the battery cells and nuts, the shape of the generated trajectories shows deviation, as shown in the purple trajectories in Figure 19a. To address the issue, the Bi-DMP previously used for unscrewing nuts needs to be updated. The 1D-SEG and G-PRM are adopted to generate new ${Tra}_{gen}$ for battery cell removal based on the yellow trajectory as shown in Figure 17b. The red trajectory in Figure 19c presents the new ${Tra}_{gen}$-based complete demonstration for removing the battery cell, and the new target point for this demonstration is positioned at the geometric centre of B6, as illustrated in Figure 19b. The light blue trajectories are generated by Bi-DMP for battery cell removal.

Figure 20 shows the implementation of the battery cell removal process. In summary, the above case studies address two critical tasks in EV battery disassembly. It can be concluded that the algorithms proposed in this paper can meet the requirements of these applications.

5. Discussion

5.1. Time Complexity Analysis

Time complexity can be used to measure the complexity of an algorithm and serves as a reference when computational resources are limited or real-time performance is required. Compared to trajectory segmentation based on neural networks, the discrete curvature-based method offers a faster computation speed because it does not require the collection of large amounts of data for training. At the same time, the Gaussian distribution, as the most widespread form of data distribution in nature, is used in the roadmap construction process in PRM. The fundamental algorithms involved in this paper are DMP, PRM, and BLS. Gaussian mixture model (GMM), TP-GMM, and RRT have similar functionalities. Analysis of trajectory segmentation algorithms using a CNN as an example [29]. The $o(·)$ is used to denote the time complexity function. According to the studies on time complexity by Liang et al. [37], Curry et al. [38], and Bianchini et al. [39], the time complexity of the aforementioned algorithms is qualitatively expressed in

Table 3. Time complexity analysis for different algorithms.

Algorithms	Time Complexity	Notes
CNN-based	$o(\mathrm{N}·\mathrm{D}·{\mathrm{F}}^2·\mathrm{C}·\mathrm{L})$	$\mathrm{N}$: the number of samples, $\mathrm{D}$: the dimension of input data (feature number), $\mathrm{F}$: the size of convolution kernel, $\mathrm{C}$: the number of convolution kernels, $\mathrm{L}$: the number of layers, $\mathrm{I}$: the number of iterations of expectation–maximisation (EM) algorithm, $\mathrm{K}$: the number of Gaussian clusters, $\mathrm{k}$: the number of sampling knots, $\mathrm{n}$: the number of sampling points, $\mathrm{E}$: the number of BLS nodes
TP-GMM	$o(\mathrm{N}·\mathrm{K}·\mathrm{D}·\mathrm{I})$
GMM	$o(\mathrm{N}·\mathrm{K}·\mathrm{D}·\mathrm{I})$
BLS	$o(\mathrm{N}·\mathrm{D}·\mathrm{E})$
DMP	$o(\mathrm{N}·\mathrm{D})$
RRT	$o(\mathrm{k}·\mathrm{l}\mathrm{o}\mathrm{g}(\mathrm{k}\left)\right)$
PRM	$o(\mathrm{n}·\mathrm{l}\mathrm{o}\mathrm{g}(\mathrm{n}\left)\right)$

.

The sorting in Table 3 is based on the time complexity of the algorithms from high to low, with a detailed analysis as follows: the convolution operations and the large number of parameters make the time complexity of the CNN the highest, especially in large-scale datasets and deep learning structures. Table 3 lists the time complexity of only a one-layer CNN. In actual computations, the time complexity of all CNN layers in the network must be summed. The time complexity of GMM is influenced by the number of Gaussian clusters [1], the total number and dimensionality of the samples, and the number of iterations of the EM algorithm. When TP is added, the data dimensionality $\mathrm{D}$ will increase, making the complexity of TP-GMM slightly higher than that of GMM with the same parameters. The time complexity of the BLS mainly depends on the number of samples, the data dimensionality, and the number of nodes in the BLS. Although the computation of the BLS relies on large-scale matrix operations, its time complexity is lower than that of a CNN. The time complexity of DMP is mainly determined by the number of samples and data dimensionality. The time complexities of PRM and RRT are related to the number of sampling points/knots. Their time complexities are similar, while PRM has an advantage in trajectory generation due to the construction of the roadmap.

5.2. Trajectory Generation Analysis among Different Algorithms

The most prominent LfD methods are the GMM-based and DMP-based algorithms.

Table 4. Quantitative comparison of various algorithms.

Algorithms	Success Rate	Running Time	Maximum Error (mm)	Difference Range (mm)
GMM	7.14%	0.97s	11.21	132.15–180.41
TP-GMM	14.29%	1.36s	105.34	201.96–970.53
DMP	57.14%	0.16s	2.76	128.65–147.88
Bi-DMP	100%	0.67s	1.48	133.74–165.79

discusses the differences among the four algorithms: GMM, TP-GMM, DMP, and Bi-DMP, during the execution of the unscrewing nuts task depicted in Figure 16. In model construction, the GMM-based algorithm requires a larger number of demonstrations, whereas DMP requires only one demonstration. Therefore, in terms of cost, DMP is the more optimal choice.

In the task of unscrewing nuts, there are differences in the performance of GMM and TP-GMM compared with DMP and Bi-DMP. GMM can generalise the distribution of demonstrations and then subsequently generate trajectories via Gaussian mixture regression (GMR) to perform tasks. However, since GMM cannot adapt to dynamic environments, it can only perform the demonstrated tasks. It means that it can only handle unscrewing a nut from a fixed position. The TP-GMM gains generalisation capabilities by adding frames of reference (frames). However, due to the sensitivity of the frames’ orientation and the inherent uncertainty of the probabilistic model, the success rate of TP-GMM is sub-optimal in high-precision applications. In addition, due to generalisation issues, DMP cannot complete all the tasks of unscrewing nuts. In contrast, DMP enhanced with a BLS can accomplish all the tasks, achieving the highest success rate. Due to differences in programming logic and computer hardware configuration, algorithm running times vary. Overall, GMM-based methods have longer running times than DMP-based methods. This is because the K-means algorithm and the EM algorithm used in GMM computations require multiple iterative operations, especially the EM algorithm. As previously mentioned, GMM and TP-GMM are probabilistic models. Thus, their generated trajectories have uncertainties, with maximum errors and difference fluctuations greater than those of DMP. The difference is the sum of the coordinates of the Euclidean distance between the generated trajectory and the demonstration trajectory for the same target. Particularly for TP-GMM, significant deviations can occur (maximum error and difference will be huge) if the frames are incorrectly specified. Therefore, for the issues discussed in this paper, DMP-based methods are more suitable.

Table 5. Quantitative comparison of PRM-based algorithms.

Algorithms	Changes in Similarity (%)
G-PRM	79.61%	88.41%	92.00%	93.39%
PRM	59.35%	No change	No change	No change

shows the differences in trajectory generation between the G-PRM and the classic PRM. As shown in the figures within the table, the blue curve represents the imitated trajectory, the purple dashed line is the trajectory generated by the classic PRM, and the red curve is the trajectory generated by G-PRM. The cyan dots are Gaussian distribution sampling points used to construct the roadmap. The blue dot is the starting point, while the red and purple dots represent the target point and the new target point, respectively. Since the classical PRM lacks imitation capabilities, the generated trajectory is merely a curve connecting the starting point and the new target point, with a similarity to the blue curve of less than 60%. Despite resampling in each iteration, the similarity to the blue curve does not improve. The G-PRM designed in this paper optimises the sampling points to continuously improve the similarity of the generated trajectory (red curve) to the blue trajectory. Its final similarity is 34.04% higher than the trajectory generated by the classical PRM.

The study has some limitations. The incremental learning process in Bi-DMP currently handles one task at a time, which may limit its efficiency in more complex applications requiring simultaneous task learning. Expanding the BLS structure to support multi-task incremental learning could significantly enhance its utility. In addition, the implementation of the method in more complex industrial applications, such as turbine and engine assembly, should be explored to fully understand its potential and limitations.

6. Conclusions

In summary, this paper proposed a systematic framework aimed at improving the application scope of DMP and BLSs. The research is organised into three main contents: (1) segmentation and extended demonstration: as a series of demonstration pre-processing methods, 1D-SEG and G-PRM are designed to segment the original demonstration and to generate ${Tra}_{gen}$. These methods effectively preserve the geometrical characteristics of the demonstration and simplify the calculation. (2) DMP modelling and incremental learning updating: Bi-DMP is trained based on ${Tra}_{gen}$ along with the force item of constructed DMP. The incremental learning mechanism of Bi-DMP enables the DMP to continually update without training from scratch. (3) EV battery disassembly case study: this study established a decommissioned battery disassembly experimental platform. Unscrewing nuts and battery cell removal have shown that the proposed algorithms effectively support industrial applications.