5.2.2.1 Related Work

HMM has been one of the most widely used methods for modeling and analysis of human motions, due to its robustness to spatiotemporal variations within sequential data (Pook and Ballard, 1993; Tso and Liu, 1997, Yang et al., 1997). On the other hand, the problem of reproduction of human motions using HMMs for task representation has been less studied in the literature. Therefore, the research presented in the book concentrates on the generation of a generalized trajectory for reproduction of HMM modeled trajectories (Vakanski et al., 2012).

Among the works in the literature that tackle the trajectory reconstruction problem, the approach reported by Tso and Liu (1997) used the observation sequence with the maximum likelihood of being generated by the trained HMM to be reproduced by the robot. However, because of the natural noise content of human movements, the demonstrated trajectories are generally not suitable for execution by the robot learner. In the work of Inamura et al. (2006), a large number of hidden state sequences for a learned HMM were synthetically generated, and a generalized trajectory was obtained by averaging over the set of trajectories that corresponded to the output distributions of the model for the given state sequences. A similar averaging approach was also employed in the article by Lee and Nakamura (2006). The drawback of the averaging approaches is that the smoothness of the reconstructed generalized trajectory cannot be guaranteed.

Another body of work proposed to generate a trajectory for reproduction of HMM modeled skills based on interpolation of trajectory key points (Inamura et al., 2003). For example, Calinon and Billard (2004) first used the forward algorithm to find the observation sequence with the maximum likelihood of the trained HMM, and then the most probable sequence of hidden states for that observation sequence was obtained by applying the Viterbi algorithm. To retrieve a task reproduction trajectory, the extracted key points were interpolated by a third‐order spline for the Cartesian trajectories and via a cosine fit for the joint angle trajectories. Although this approach generated a smooth trajectory suitable for robot reproduction, the resulting trajectory corresponded to only one observed demonstration. What is preferred instead is a generalized trajectory that will encapsulate the observed information from the entire set of demonstrations. The HMM‐based key points approach for trajectory learning was also employed by Asfour et al. (2006). The authors introduced the term of common key points (which refers to the key points that are found in all demonstrations). Thus, the generalized trajectory included the relevant features shared by the entire set of individual demonstrations. Their approach was implemented for emulation of generalized human motions by a dual‐arm robot with the generalized trajectory obtained by linear interpolation through the common key points.

The approach presented in the book is similar to the works of Calinon and Billard (2004) and Asfour et al. (2006), in using key points for describing the important parts of the trajectories for reproduction. Unlike these two approaches, a clustering technique is used for initial identification of candidate key points prior to the HMM modeling. The main advantage of the presented approach is in employing the key points from all demonstrations for reconstruction of a trajectory for task reproduction. Contrary to the work of Asfour et al. (2006), the key points that are not repeated in all demonstrations are not eliminated from the generalization procedure. Instead, by introducing the concept of null key points, all the key points from the entire set of demonstrations are considered for reproduction. Furthermore, by assigning low weights for reproduction to the key points that are not found in most of the demonstrations, only the most consisted features across the set of demonstrations will be reproduced. A detailed explanation of the approach is given in the following sections.

5.2.2.2 Modeling

5.2.2.2.2 HMM Initialization

The efficiency of learning with HMMs in a general case depends directly on the number of available observations. Since in robot PbD the number of demonstrations is preferred to be kept low, a proper initialization of the HMM parameters is very important. Here, the initialization is based on the key points assignment with the LBG algorithm. The minimum distortion criterion (Linde et al., 1980; Tso and Liu, 1997) is employed to select a demonstration X_σ, based on calculating the time‐normalized sum‐of‐squares deviation, that is,

(5.8)

where α_m(t_k) is the combined position/velocity vector of time frame k in the demonstration m, and is the centroid corresponding to the label (o_m,k).

Accordingly, the number of hidden states N_s of HMM is set equal to the number of regions between the key points of the sequence X_σ (i.e., the number of key points of X_σ plus one). The HMM configuration known as a Bakis left‐right topology (Rabiner, 1989) is employed to encode the demonstrated trajectories. The self‐transition and forward‐transition probabilities in the state‐transition matrix are initialized as follows:

(5.9)

where τ_i,σ is the duration in time steps of the state i in demonstration X_σ, and W is a stochastic normalizing constant ensuring that . The transitions allow states to be skipped in the case when key points are missing from some of the observed sequences. All other states transition probabilities are assigned to zeros, making transitions to those states impossible. The output probabilities are assigned according to

(5.10)

where n_i,σ(c_f) represents the number of times the symbol c_f is observed in state i for the observation sequence . The initial state probabilities are set equal to .

5.2.2.3 Generalization

5.2.2.3.1 Temporal Normalization of the Demonstrations via Dynamic Time Warping

Since the demonstrated trajectories differ in lengths and velocities, the extracted sequences of trajectory key points correspond to different temporal distributions across the trajectories. For generalization purposes, the set of key points needs to be aligned along a common time vector. This will ensure tight parametric curve fitting in generating a generalized trajectory through interpolation of the extracted key points. Otherwise, the temporal differences of the key points across the demonstrations could cause spatial distortions in the retrieved reproduction trajectory. The temporal reordering of the key points is performed with the dynamic time warping (DTW) algorithm through the alignment of the entire set of demonstrated trajectories with a reference trajectory.

To select a reference sequence for time warping of the demonstrations, the forward algorithm (Section 4.2.1) is used for finding the log‐likelihood scores of the demonstrations with respect to the trained HMM, that is, . The observation with the highest log likelihood is considered as the most consistent with the given model; hence, the respective trajectory X_ϕ is selected as a reference sequence.

Next, the time vector of the trajectory X_ϕ, that is, , is modified in order to represent the generalized time flow of all demonstrations. This is achieved by calculating the average durations of the hidden states from all demonstrations , for . Since the key points are taken to be the transitions between the hidden states in the HMM, the new key point time stamps for the demonstration are reassigned as follows:

(5.11)

where denotes the new time index of the jth key point, and since the first key point is always taken at the starting location of each trajectory. Subsequently, a linear interpolation between the new time indexes of the key points is applied to find a warped time vector . The reference sequence is obtained from the most consistent demonstrated trajectory X_ϕ with time vector .

Afterwards, multidimensional DTWs are run on all demonstrated trajectories against the time‐warped data sequence . It resulted in a set of warped time waves, G₁, G₂, …, G_M, that aligned every trajectory with the trajectory. Accordingly, all key points are shifted in time by reassigning the time index of every key point κ_j,m with the time index contained in the jth member of G_M. The resulting warped sequences of key points are denoted by , for .

With the key points from all demonstrations temporally aligned, parametric curve fitting across each dimension would suffice to create a generalized trajectory that represents the set of demonstrated trajectories. However, this approach would fail to separate the deliberate and precise motions from the transitional portions of the trajectories. Therefore, weighting coefficients are assigned for the different portions of trajectories, as bias coefficients of their relative importance for reproduction.

5.2.2.3.2 Key Points Weighting

Introducing the weighting coefficients was inspired by the works of Billard et al. (2004) and Calinon and Billard (2007), where the relative importance of each dimension from a set of measured variables was compared, and for reproduction purposes greater weights were imposed on the dimensions that were less variant across the different demonstrations.

Similarly, in this work weighting coefficients are assigned to all clusters of key points. The weighting is based on the assumption that the key points that are not closely clustered in comparison to the rest of the key points are generally less significant for reproduction. The weighting concept is illustrated in Figure 5.2.

The root mean‐square (RMS) error is employed to measure the closeness of each cluster of key points as follows:

(5.12)

where denotes the centroid of the cluster Q_j of key points with time index j. For each cluster of key points Q_j, a normalized weighting coefficient is assigned based on

(5.13)

where ε_min and ε_max are empirically chosen threshold RMS error values. The introduced weighting scheme in (5.13) sets the weights to 0 and 1, when the RMS error is below or above a predefined threshold value, respectively. The upper threshold value ε_max is typically chosen at two or three standard deviations of the data points from the centroids of the corresponding cluster. The lower threshold value ε_min is typically set equal to a half of the standard deviation of the data points, in order to obtain high weights for reproduction of the data points that are closely clustered. The selection of the threshold values ε_min and ε_max for the particular trajectories used for evaluation of the presented method is explained in Section 5.2.2.4.

It should be noted that the null key points are not used for weighting or interpolation. Instead, they were introduced to ensure that the key points with the same index correspond to identical spatial positions across the trajectories.

5.2.2.4 Experiments

A block diagram that describes the high‐level information flow for the developed approach is presented in Figure 5.3. The algorithm automatically produces a generalized version of the set of demonstrated trajectories, based on a series of data processing steps. A human operator intervened in collecting the demonstrated trajectories, and in the transfer of the derived strategy for task reproduction to the robotic platform for execution.

To evaluate the approach, two experiments are conducted for a manufacturing process of painting. A hand tool is moved by an operator in a simulated painting process, by traversing the areas of the parts that are required to be painted. The first experiment follows a simple trajectory, while the second experiment entails a more complex geometry and trajectory that includes nonplanar regions and waving motions of the painting tool with different amplitudes.

The experiments were conducted in the National Research Council (NRC)—Aerospace Manufacturing Technology Centre (AMTC) in Montréal, Canada.

5.2.2.4.1 Experiment 1

The experimental setup for the first experiment is shown in Figure 5.4. A demonstrator is performing virtual painting of a panel using a hand tool that serves as a spray gun. The attached optical markers on the tool are tracked by the optical tracking system Optotrak (Section 2.1) shown in Figure 5.4b. The task trajectory consists of (i) approaching the upper left‐side corner of the panel from an initial position of the tool, (ii) painting the rectangular contour of the panel in a clockwise direction, (iii) painting the inner part of the panel with waving left‐to‐right motions, and (iv) bringing the tool back to the initial position.

The task is demonstrated four times by four different operators, that is, the measurement set consists of trajectories X_m, for , containing the position and orientation data of the tool with respect to a reference frame (Figure 5.4c). The robustness of the learning process can be enhanced by eliminating the demonstrations that are inconsistent with the acquired dataset (Aleotti and Caselli, 2006). Otherwise, the learning system might fail to generalize correctly from a set of demonstrations which are too dissimilar. The sources of variations can be attributed to the uncertainties arising from (i) involuntary motions of the demonstrator, (ii) different conditions in performing the demonstrations, (iii) difference in the performance among the demonstrators, (iv) fatigue due to multiple demonstrations of the same task, etc. The presented work emphasizes on the problem of inconsistency within a set of demonstrations executed by different human demonstrators, through the introduction of a preliminary step for evaluation of the demonstrators’ performance.

This step carries out an initial comparison of the distributions of the first, second, and third derivative of position with respect to time, that is, velocities, accelerations, and jerks, respectively, of the demonstrated trajectories. For each subject, the values of the parameters from all trajectories are combined, and the distributions of the parameters are represented with box plots in Figure 5.5. The plots, for example, indicate that the demonstrations by Subject 4 have the greatest deviations of the accelerations and the jerks, while for the rest of the subjects the distributions are more uniform. Therefore, the demonstrations by Subject 4 were eliminated from the recorded dataset. These parameters can also be exploited as criteria for evaluation of the demonstrations with regard to the motor capabilities of the robot learner. For instance, if the levels of accelerations of the demonstrations are greater than the maximum accelerations available from the robot links motors, the executed task will differ from the demonstrated task, and it can even lead to a task failure, depending on the circumstances. Consequently, the input set of trajectories is reduced to trajectories, demonstrated by Subjects 1, 2, and 3. This initial refinement results in a more consistent initial dataset, and it leads to generation of a more robust task model. The Cartesian positions of the demonstrations are shown superimposed in Figure 5.4c. The number of measurements of the demonstrations (i.e., T_m) varies from 3728 to 7906 measurements, with an average of 5519. One can infer that despite the similarity of the trajectories’ geometry in Figure 5.4c, the velocities of the demonstrations varied significantly between the subjects.

After an initial preprocessing of the data, which included smoothing of the recorded trajectories with a Gaussian filter of length 10, candidate key points are identified using the LBG algorithm, as explained in Section 5.2.1. The parameter ς, related to the delay of time frames for the velocity calculation, is chosen equal to 20 sampling periods (i.e., 200 milliseconds), and the number of clusters ρ is set to 64. These parameters are chosen heuristically to balance the complexity and the duration of the demonstrated trajectories. The number of key points per trajectory ranges from 77 to 111 (average value of 87.75 and standard deviation of 11.7). For the minimum distortion trajectory X₁₂ (found using (5.8)), 82 key points are identified. The key points which are closer than 20 time frames are eliminated from the minimum distortion trajectory, resulting in 72 key points in total (shown in Figure 5.6).

Afterward, the recorded trajectories are mapped into discrete sequences of observed symbols, and discrete HMMs are trained. For that purpose, Kevin Murphy’s HMM toolbox (Murphy, 1998) is used in MATLAB environment. The number of discrete observation symbols N_o is set to 256. The number of hidden states N_s is set equal to the number of key points of the minimum distortion trajectory plus 1, that is, 73. The Viterbi algorithm is employed for finding the sequences of hidden states for each observation sequence, and key points are assigned to the position and orientation values of the trajectories that correspond to the transitions between the hidden states in the observation sequences. Thus, including the starting and ending points of the trajectories, each demonstration was described by a time‐indexed sequence of 74 key points. As explained earlier, for the observation sequences with missing hidden states, the coordinates of the corresponding missing key points are set to zero (for that reason, they were referred to as null key points). This procedure results in spatiotemporal ordering of the key points for all trajectories, with a one‐to‐one mapping existing between the key points across the different demonstrations (e.g., for each trajectory, the key point with the index corresponds to the upper right corner of the contour). Log likelihoods of the observation sequences for the trained HMM are obtained with the forward algorithm. The observation sequence had the highest log likelihood; therefore, the demonstration X₃ is used for temporal alignment of the set of trajectories with the DTW algorithm.

To obtain a generalized trajectory that would represent the entire set, its length should be approximately equal to the average length of the trajectories in the set. Hence, the average durations of the hidden states for the trained HMM are calculated, whereby the newly created sequence of hidden states has length of 5522 time frames. This value is indeed very close to the average length of the demonstrations, being 5519. This time vector is used to modify the time flow of the demonstration X₃ as explained in Section 5.2.2.3.1. Subsequently, the resulting sequence is used as a reference sequence for the DTW algorithm, that is, to align each of the demonstrated trajectories against it. The DTW alignment resulted in time‐warped trajectories with length of 5522 time frames, from which the time instances of the key points are extracted. The Cartesian positions coordinates x, y, and z of the DTW‐reassigned key points from all demonstrations are shown in Figure 5.7.

For the approaching and departing parts of trajectories in Figure 5.7, the variance among the key points is high, which implies that these parts are not very important for reproduction. The distributions of the RMS errors for the clusters of key points along x, y, and z coordinates are shown in Figure 5.8a. There are 73 clusters of key points in total, where the starting and ending subsets of key points, which correspond to the approaching and departing parts of the trajectories, exhibit larger deviations from the centroids of the clusters. Therefore, weighting coefficients with lower values are assigned for those parts of the trajectories. The variance of the key points associated with painting the contour and the inner side of the panel is low, and thus the computed weighting coefficients have values close to 1. Values of 0.5 and 2 standard deviations of the RMS errors for the corresponding cluster are adopted for ε_min and ε_max, respectively, in (5.13). The weights for such adopted thresholds are displayed in Figure 5.8b. For comparison, the weighting coefficients for the key point clusters with the threshold values of 1/6 and 6 standard deviations are shown in Figure 5.8c. In this case, the weights distribution will result in a reproduction trajectory that will fit more closely the demonstrated trajectories. However, such choice of the parameters will not have very significant impact on the final outcome, that is, the approach is not highly sensitive to the choice of the threshold values ε_min and ε_max.

The last step entails smoothing spline interpolation through the key points by using an empirically determined smoothing factor , which gives both satisfactory smoothness and precision in position for the considered application. For the orientation of the tool expressed in Euler roll–pitch–yaw angles, the 1D generalization results are given in Figure 5.9. The plots show the angle values of the key points and the interpolated curves for each of the angle coordinates. A plot of the Cartesian coordinates of the resulting generalized trajectory X_gen is shown in Figure 5.10.

For comparison, the distributions of velocities, accelerations, and jerks for the demonstrated trajectories by the three subjects and the generalized trajectory are shown in Figure 5.11. First, from the plot of velocity distributions, it can be concluded that the mean value of the generalized trajectory is approximately equal to the average value of the velocities from the demonstrations by all subjects. This is the result of the generalization procedure, that is, the demonstrations by all subjects are used to find a generalized trajectory that will represent the entire set of demonstrations. Second, the distributions of accelerations indicate that the values of the generalized trajectory are lower than the values of the demonstrated trajectories. It is due to the fact that spline smoothing interpolation is employed for reconstruction of the generalized trajectory, which reduces the changes of velocities present in human motions. Third, the jerk distribution for the generalized trajectory is significantly condensed when compared to the demonstrated trajectories. Namely, the inconsistency of the human motions represented by high levels of jerk is filtered out by the generalization procedure, resulting in a smooth trajectory with low levels of jerk.

5.2.2.4.2 Experiment 2

The second experiment entails painting a panel with more complex geometry, shown in Figure 5.12a. The objective is to paint the top side of the panel, and the horizontal plane on the right‐hand side of the panel. Both areas are shown bordered with solid lines in Figure 5.12a. The three optical markers shown in the picture are used to define the reference coordinate system of the panel. The task is demonstrated five times by a single demonstrator, and all trajectories are shown superimposed in Figure 5.12b. Different from Experiment 1 where the shape of the generalized trajectory for the painting part is obvious from Figure 5.4c, in the case of Experiment 2 it is not easy to infer the shape of the generalized trajectory, especially for the small right‐hand plane where the waving motions are not controlled.

The lengths of the demonstrations range between 4028 and 4266 measurements. The initial selection of candidate key points with the number of clusters ρ set to 64 results in numbers of key points between 122 and 134 per trajectory. Based on the sequence with minimum distortion, the number of hidden states N_s for training the HMM is set to 123. Afterward, the demonstrations are temporally aligned with DTW, with the length of the sequence of average states duration equal to 4141 time frames. The resulting set of key points from all five demonstrations is interpolated to generate a generalized trajectory X_gen (Figure 5.13a). The plot also shows the approach and retreat of the painting tool with respect to the panel.

The generalized trajectory for the pose of the painting tool was transferred to a desktop robot (Motoman^® UPJ), and afterward it was executed by the robot (Figure 5.13b). The height of the robot’s end‐effector was set up to 30 millimeters above the part. The objective of the experiments was not to perform the actual painting with the robot, but to illustrate the approach for learning complex trajectories from observations of demonstrated tasks. The results demonstrate that with encoding of demonstrations at the trajectory level, complex demonstrated movements can be modeled and generalized for reproduction by a robot in a PbD environment.

5.2.2.5 Comparison with Related Work

To compare the presented method with the similar works in the literature, the obtained trajectory for task reproduction based on the approach reported by Calinon and Billard (2004) is shown in Figure 5.14. The demonstration set for Experiment 1 from the Section 5.2.2.4.1 is used, and the result is obtained by using the same parameter values reported in the aforementioned section. The initial key points selection based on LBG algorithm is employed, in order to compare the approaches on the most similar basis. The trajectory for task reproduction is generated by smoothing spline interpolation of the key points generated from the Viterbi path of the most consistent trajectory, as explained in the work of Calinon and Billard (2004). The most consistent trajectory is shown with the dashed line in Figure 5.14. The main drawback of this approach is that it does not generalize from the demonstrated set. Instead, the reproduction trajectory is obtained by using a criterion for selecting the best demonstration from the recorded set, followed by a smoothing procedure. On the other hand, in the approach presented in the book the trajectory for the task reproduction (shown in Figure 5.10) is obtained by taking into account important features from the entire set of demonstrated trajectories, and it does not rely so heavily on only one of the demonstrated trajectories.

Another advantage of the presented approach with respect to the work by Calinon and Billard (2004) is in providing a mechanism for interpretation of the intent of the demonstrator. Based on the variance of the key points corresponding to the same spatial features across the demonstrated set, the method presented here assigns different weights for reproduction purposes. For instance, for the approaching and departing (transitional) sections of the trajectories, the variability across the demonstrated set is high, and subsequently the generalized trajectory does not follow strictly the trajectories key points. On the other hand, the section of the trajectories corresponding to painting the panel is more constrained, which is reflected in the generation of the generalized trajectory with tight interpolation of the identified key points. It can be observed in Figure 5.10 that the approaching and departing sections are smoother, due to the higher variance across the demonstrations, whereas in Figure 5.14 the approaching and departing sections simply follow the most consistent trajectory. The procedure adopted in the book has similarities to the work of Calinon (2009) based on GMM/GMR, where the sections of the demonstration with high values of the Gaussian covariances are considered less important for reproduction purposes, and vice versa.

The approach of Asfour et al. (2006), which employs generalization from the set of demonstrations by using the concept of common key points, is similar to the presented approach, although it may provide suboptimal results in some cases. For instance, if a key point is missing in one of the demonstrated trajectories, and it is present in all of the remaining trajectories, this key point will be eliminated as not being common and will not be considered for generalization purposes. To avoid this case, the approach might require applying a heuristic for analysis of the demonstrations, and elimination of the suboptimal demonstrations.

The results from the developed method and the methods reported in Calinon and Billard (2004) and Asfour et al. (2006) are compared by using the RMS error as a metric for evaluation of learned trajectories (Calinon et al., 2010)

(5.14)

where m₁ and m₂ are used for indexing the trajectories. The trajectories were scaled to sequences with the same length by using two techniques: linear scaling and DTW alignment. With the linear scaling technique, the length is set equal to the average length of the demonstrated trajectories. For the DTW scaling, the entire set of trajectories is aligned against one of the demonstrated trajectories. Rather than selecting the generalized trajectory for scaling the dataset with the DTW, some of the demonstrated trajectories were used as reference sequences for alignment, to avoid biased results for the RMS differences.

A color map chart with the RMS differences for the linearly scaled trajectories for Experiment 1 is shown in Figure 5.15. The first three trajectories in the chart represent the generalized trajectories obtained by the presented approach—X_G1, by the approach proposed in Calinon and Billard (2004)—X_G2, and the one proposed in Asfour et al. (2006)—X_G3. The rest of the cells correspond to the demonstrated trajectories X₁–X₁₂. The darker colors in the chart pertain to cells with smaller RMS differences. For instance, the trajectory X₅ is the most dissimilar to the whole set, and it can be noted that the trajectories demonstrated by the same subject X₁–X₄, X₅–X₈, and X₉–X₁₂ are more similar when compared to the trajectories demonstrated by the other subjects.

The sums of the RMS differences for the trajectories from Figure 5.15 are given in Figure 5.16a. The RMS differences for scaling of the trajectories with DTW are also displayed. The demonstrated trajectories X₄, X₅, and X₁₂ are chosen as reference sequences for aligning the trajectories. The DTW algorithm is repeated for three of the demonstrated trajectories to avoid any subjective results that can arise from using a particular trajectory as a reference sequence. As expected, the DTW algorithm aligns the trajectories better than the linear scaling technique. For the presented four cases, the cumulative RMS differences for the trajectory obtained by the presented approach (X_G1) are the smallest, which infers that it represents the best match for the demonstrated set.

Similarly, for Experiment 2 from the Section 5.2.2.4.2, the sums of the RMS differences are presented in Figure 5.16b. The DTW alignment is repeated for two of the demonstrated trajectories, X₂ and X₄. Based on the presented results, it can be concluded that the generalized trajectory generated by the described approach fits better the demonstrated set, when compared to the similar approaches proposed by Calinon and Billard (2004) and Asfour et al. (2006).

Another work in the literature that has some similarities to the one presented in the book has been reported by Gribovskaya and Billard (2008). In their study, DTW was employed initially to align the recorded sequences. Afterward, a continuous HMM was employed for encoding the trajectories. One drawback of such procedure is that the DTW algorithm distorts the trajectories. As a result, important spatial information from the demonstrations may be lost. In addition, the velocity profile of the DTW aligned demonstrations is deformed. In the method presented here, HMM is used for segmentation of the trajectories as recorded, and DTW is used before the interpolation procedure, only to shift the key points to a trajectory with a given time duration (which is found as the average of time durations of the hidden states form HMM). Hence, the velocity information is preserved for modeling the demonstrations with HMM, and the DTW affects the trajectories only at the level of key points.

The computational expense of the presented approach is compared with the state‐of‐the‐art approach reported in Calinon (2009), which uses GMM and GMR for modeling and generalization of continuous movements in a robot PbD setting. The simulations are run on a 2.1 GHz dual‐core CPU with 4 GB of RAM running under Windows XP in the MATLAB environment. Since there are slight variations for the processing times in different simulation trials, the codes were run five times, and the mean values and standard variations are reported. The computation times for obtaining a generalized trajectory from the sets of measured trajectories for both experiments from Section 5.2.2.4 are presented in Table 5.2. The DTW alignment step of the method is performed with a MATLAB executable (MEX) file, which increases the computational speed of the DTW algorithm for about 6–8 times compared to the standard MATLAB m‐files. The results of the computational cost for the presented approach displayed in Table 5.2 indicate that most of the processing time is spent on the learning of HMM parameters and on the discretization of the continuous trajectories.

The computation times for the GMM/GMR method are reported in Table 5.3, with using the same MATLAB MEX subroutine for performing the DTW alignment phase. For Experiment 1, the trajectory X₁₂ (with the length of 4939 measurements) is selected as a reference sequence for aligning the set of 12 demonstrated trajectories, whereas for Experiment 2 the trajectory X₂ (4048 measurements) is selected as a reference sequence for aligning the set of five trajectories. For the steps of GMM encoding and GMR generalization, the MATLAB codes provided in Calinon (2009) are used. The number of Gaussian components for modeling the trajectories was approximated by trial and error, and the values used for the reported times are 40 for the first experiment and 75 for the second experiment.

The generalized trajectories generated from the GMM/GMR method are given in Figure 5.17. Based on the results from Table 5.3, it can be concluded that the computational cost of modeling the 6D input data into a large number of Gaussians is computationally expensive. The total computation times with the presented approach are less than the times obtained with the GMM/GMR approach for both sets of trajectories.

5.2.3.1 Related Work

One of the application domains where CRFs have been most extensively utilized is the language processing. Examples include part‐of‐speech tagging (Lafferty et al., 2001; Sutton and McCallum, 2006), shallow parsing (Sha and Pereira, 2003), and named‐entity recognition (McDonald and Pereira, 2005; Klinger et al., 2007). CRFs were reported to outperform HMMs for classification tasks in these studies. Other areas of implementation include image segmentation (He et al., 2004; Quattoni et al., 2005), gene prediction (DeCaprio et al., 2007), activity recognition (Sminchisescu et al., 2005; Vail et al., 2007a), generation of objects’ trajectories from video sequences (Ross et al., 2006), etc.

In the robot PbD environment, Kjellstrom et al. (2008) presented an approach for grasp recognition from images using CRF. Martinez and Kragic (2008) employed the support vector machine (SVM) algorithm for activity recognition in robot PbD by modeling the demonstrated tasks at a symbolic level, followed by using CRF for temporal classification of each subtask into one of the several predefined action primitives classes. The comparison results of recognition rates between HMM and CRF indicated similar performance of both methods for short actions, and higher recognition rates by CRF for the case of prolonged continuous actions.

CRF has not been implemented before for acquisition of skills at the trajectory level in the area of robot PbD. Therefore, Section 5.2.3.2 presents one such approach (Vakanski et al., 2010). Similarly to the approach presented in Section 5.2.2 that employed HMM for encoding trajectories, the hidden states of the demonstrator’s intention in executing a task are related to the trajectory key points in the observed sequences. CRF is employed for conditional modeling of the sequences of hidden states based on the evidence from observed demonstrations. A generalized trajectory for reproducing the required skill by the robot is generated by cubic spline interpolation among the identified key points from the entire set of demonstrated trajectories.

5.2.3.2 Feature Functions Formation

The developed framework for robot learning from multiple demonstrations employs the optical motion capture device Optotrak Certus^® (NDI, Waterloo, Canada) for the task perception (Section 2.1). Hence, the observation sequences are associated with the continuous trajectories captured during the demonstration phase. This section discusses the formation of the feature functions in a linear chain CRF structure (Section 4.3.1), for the case when the observation sequences are human demonstrated trajectories. In general, CRFs can handle observed features of either continuous or discrete nature, whereas the hidden state variables are always of discrete nature. Structuring of the feature functions for both continuous and discrete observations is discussed next.

A CRF model designed for classification of human activities from continuous trajectories as observation data is presented in the work of Vail et al. (2007b). The authors introduced the following continuous‐state feature functions:

(5.15)

for , where N_s denotes the number of hidden states, and d is used for indexing the dimensionality of the observed Cartesian positions o_k. The notation pertains to a binary indicator function, which equals to 1 when the state at time t_k is i, and 0 otherwise. To enhance the classification rate of the CRF in the article of Vail et al. (2007b), additional continuous‐state feature functions are added to the model, corresponding to the observed velocities, squared positions, etc.

The transition feature functions for encoding the transition scores from state i at time to state j at time t_k are defined as follows:

(5.16)

for .

An important remark reported in the article (Vail et al., 2007b) is that the empirical sum of the features in the log‐likelihood function (4.23) can be poorly scaled when continuous features functions are used. This can cause slow convergence of the optimization algorithms, or in the worst case, the algorithm would not converge. As a remedy, the authors proposed to perform normalization of the continuous observations to sequences with zero mean and unity variance, which improved the classification results.

Nevertheless, most of CRF’s applications in the literature deal with categorical input features (e.g., labeling words in a text), rather than continuous measurements. In that case, the observed features are mapped to a finite set of discrete symbols . The state feature functions are defined as binary indicators functions for each state‐observation pairs as follows:

(5.17)

for and .

The transition feature functions are defined in an identical manner as in (5.16). The structure of CRFs allows additional observed features to be easily added to the model by generalizing the feature functions. For instance, the dependence of state transitions to the observed symbol can be modeled by adding additional feature functions .

5.2.3.3 Trajectories Encoding and Generalization

The implementation of CRF in this work is based on discrete observation features, due to the aforementioned scalability problems with the continuous features. For that purpose, the continuous trajectories from the demonstrations are mapped to a discrete set of symbols by using the LBG algorithm, as described in Section 5.2.1. The output of the LBG clustering procedure is a set of discrete sequences .

Training the CRF model requires provision of the corresponding sequences of hidden states (i.e., labels) for each observation sequence . Similar to the approach described in Section 5.2.2, the labels are assigned to the key points along the trajectories, which are characterized by significant changes in position and velocity. The initial labeling of the trajectories is performed with the LBG technique. The clustering is carried over a set of concatenated normalized positions and velocities vectors from all demonstrated trajectories, that is, for , , and . The transitions between the cluster labels are adopted as key points. The detected initial key points resulted in the labeled sequences denoted as .

A common challenge associated with the task modeling in robot PbD is the limited number of demonstrations available for estimation of the model parameters (since it may be frustrating for a demonstrator to perform many demonstrations, and in addition, the quality of the performance can decrease due to fatigue or other factors). Therefore, to extract the maximum information from limited training data, as well as to avoid testing on the training data, the leave‐one‐out cross‐validation technique is adopted for training purposes. This technique is based on using a single observation sequence for validation, and using the remaining observation sequences for training. For instance, to label the observed sequence which corresponds to the third demonstration , a CRF model is trained on the set consisting of the observed sequences , and the corresponding state sequences obtained from the initial selection of the candidate key points. The set of parameters Λ₃, estimated during the training phase using (4.23) and (4.24), is afterward employed to infer the sequence of states which maximizes the conditional probability . This procedure is repeated for each observed sequence from the dataset.

The inference problem of interest here consists of finding the most probable sequence of labels for an unlabeled observation sequence. The problem is typically solved either by the Viterbi method or by the method of maximal marginals (Rabiner, 1989). The maximal marginals approach is adopted in this work. This approach is further adapted for the task at hand by introducing additional constraints for the state transitions. Namely, the hidden state sequences in the presented CRF structure are defined such that they start with the state 1 and are ordered consecutively in a left‐right structure, meaning that the state 1 can either continue with a self‐transition or transition to the state 2. Furthermore, some states might not appear in all sequences, for example, a sequence can transition from the state 5 directly to the state 7, without occurrence of the state 6. By analogy with the left‐right Bakis topology in HMMs, the contributions from the potential functions which correspond to the transitions to the past states (i.e., in (4.20) for and ) are minimized by setting low values for the parameters . Thereby, the transitions to the previously visited states in the model are less likely to occur. Additionally, another constraint is introduced that minimizes the possibility of transitioning to distant future states, by setting low values for the parameters for . As a result, the potential functions corresponding to the transitions to more than two future states have approximately zero values. These constraints reflect the sequential ordering of the hidden states in the CRF model for the considered problem.

The estimation of the hidden state at time t_k is based on computing the maximal marginal probabilities of the graph edges, and are solved by the forward–backward algorithm

(5.18)

where and denote the forward and backward variables which are calculated in the same recursive manner as with HMMs (Section 4.2.1). The functions Ω_k in (5.18) correspond to the transition feature functions in (4.20) (Sutton and McCallum, 2006), that is,

(5.19)

for .

The classification rate comparing the initial sequence of states and the sequence of states obtained by the CRF labeling is used as a measure of fit for each trajectory m

(5.20)

The generation of a trajectory for reproduction of the demonstrations is analogous to the generalization procedure for HMM described in Section 5.2.2.3. To recapitulate, for each sequence of hidden states, the key points are taken to be the transitions between the states. The temporal variance across the set of sequences is tackled by employing the multidimensional DTW algorithm. The trajectory X_ϕ corresponding to the states sequence with the maximum classification rate in (5.20) is selected as a reference sequence, against which all the other trajectories are aligned. This process results in temporally warped key points sequences , where the key points’ time stamps correspond to the time vector of the reference trajectory X_ϕ. Weighting coefficients are introduced based on the variance of key points across the demonstrations associated with the same time stamps, in order to account for the relative importance for interpolation. The key points from all trajectories are afterward interpolated with a smoothing cubic spline, yielding a generalized trajectory X_gen.

5.2.3.4 Experiments

The presented concept for robot PbD is assessed through simulations of two types of manufacturing processes: painting and shot peening. Perception of the demonstrated trajectories is performed with an optical tracking device. The approach is implemented in MATLAB environment by using the CRF chain toolbox (Schmidt et al., 2008).

5.2.3.4.1 Experiment 1

The experimental setup shown in Figure 5.18a is employed to simulate the trajectories of a tool for painting target objects. The painting task is demonstrated 14 times, and it consists of first painting the contour followed by painting the interior part of the panel. The demonstrations are shown in Figure 5.18c. The trajectories lengths T_m range from 1393 to 2073 time frames. For this task, the orientation of the tool is required to be normal to the panel during the painting; therefore, it is not considered as a relevant discriminative feature. For the initial selection of key points, the trajectories are automatically segmented with the LBG algorithm using 32 clusters of discrete vectors. A total of 33 transitions between the LBG’s discrete symbols appeared in all trajectories, and these are adopted as the initial key points. The first and the last measurements of each trajectory are added as key points, resulting in 35 key points per trajectory. One sample trajectory with the identified key points is shown in Figure 5.18b. The demonstrated trajectories are discretized by utilizing the vectors consisting of normalized 3D positions and velocities . With 32 symbols used per feature, the total number of observation symbols N_o amounted to 192.

The percentage of the average classification rates of CRFs from (5.20) calculated for different observation features are given in Table 5.4. The best classification rates are obtained when both velocities and positions are used as observation features. The generalized trajectory is shown by the thick line in Figure 5.18c, superimposed with the demonstrated trajectories. Note that even if some key points were wrongfully classified by the CRF, they would not have big impact on the generalization due to the introduced weighting scheme, which assigns low weights for fitting to the key points with high variability relative to the key points from the other trajectories.

	Classification rates
1. Position	81.17 (±9.20)
2. Velocity	76.86 (±3.27)
3. Position and velocity	89.11 (±3.13)

5.2.3.4.2 Experiment 2

The second task for evaluation of the presented approach pertains to application of a shot peening procedure for surface finishing (Irish et al., 2010). Shot peening is a process of impacting material’s surface with small spherical elements called “shots.” This process produces a compressive residual stress layer on material’s surface, which increases the resistance to cracks caused by material fatigue or stress corrosion.

Demonstration of the peening task is performed seven times, with the number of measurements T_m varying between 2438 and 2655. Conversely to the painting process in Experiment 1, where the demonstrations were simulated in the university’s lab, this task’s demonstrations were captured for a process of peening in a real industrial environment. The trajectories represent waiving motions over a curved surface (Figure 5.19). Tool’s positions, velocities, and orientations are examined as observation features. The number of clusters ρ for the initial key points selection with the LBG algorithm is set to 4. This choice is appropriate for the simple waving pattern of the task. A total of 41 key points per trajectory are identified, making the number of hidden states N_s equal to 40. One sample trajectory with the initially selected key points is shown in Figure 5.19a. The discretization of the trajectories is performed with the number of pre‐specified clusters for LBG equal to 16, that is, there are 16 discrete observation symbols for each dimension of the observed sequences.

The classification rates of the trajectories segmentation with CRF are given in Table 5.5. From the results provided, it can be concluded that the velocity features are the most suitable for CRF segmentation of this example. The reason behind it is that the key points corresponding to the peaks and troughs of the traversed waving motions (Figure 5.19) are associated with the minimums of the trajectories velocities, whereas the middle parts of the waving motion corresponding to the waves’ baseline are associated with the maximums velocities. Hence, these features are invariant for the entire set, and are informative for classification purposes. With regard to the tool’s orientation, only the angles around the axis normal to the surface are taken into consideration, since the orientations around the other two axes are almost constant and do not convey important information for classification purposes. The results provided in Table 5.5 indicate that the orientation information does not improve the classification rates. Moreover, the CRF classification with the orientation as the only observation feature failed. The generalized trajectory for the case of position and velocity employed as observation features is shown in Figure 5.19b. The demonstrated trajectories are not shown in Figure 5.19b, since the demonstrations are too dissimilar, and it would be impossible to differentiate the individual demonstrations in a superimposed plot.

	Classification rates
1. Position	79.06 (±4.33)
2. Velocity	94.30 (±1.79)
3. Position and velocity	94.70 (±1.55)
4. Position, velocity, and orientation	93.49 (±4.37)

5.2.3.5 Comparisons with Related Work

The presented approach is compared to HMM, due to its wide usage for labeling and analysis of sequential data. To ensure that both methods are compared on an equal basis, a discrete form of HMM is employed, where the input data are the same discrete observation sequences and the label sequences , which were used in the CRF approach. The multidimensional HMM reported by Yang et al. (1994) is implemented. Different from the 1D HMM, the multidimensional HMM has one observation matrix B^(d) for each dimension d of the observed data. The model is built on an assumption that each dimension of the observed sequences is independent from the other dimensions. On the other hand, the state transition matrix A has the same form as in the 1D case. The training and inference problems are solved with small modifications of the 1D HMM (see the work of Yang et al. (1994) for details). To maintain consistency with the adopted cross‐validation testing for CRF, the HMM parameters are initialized using the CRF‐labeled data for one of the trajectories (e.g., and ), followed by the model parameters estimation with the Baum–Welch algorithm for the rest of the trajectories (e.g., ). The most probable sequence of hidden states for each observation trajectory and the given model are obtained by the Viterbi algorithm. This procedure is repeated for each of the observed trajectories.

For Experiment 1 related to the painting task, the mean classification rates for the multidimensional HMM obtained from (5.20) are reported in Table 5.6. For comparison, the classification rates obtained by the CRF approach from Table 5.4 are also provided in the rightmost column. For the three cases of observation features considered, the HMM produced lower number of correct labels than the CRF. The results from statistical analysis using paired t‐tests showed statistically significant difference in the classification rates by the CRF and HMM at a 0.05 level. The following p‐values were obtained for the different observation features: position (), velocity (), and position/velocity ().

For the peening task in Experiment 2, the classification rates obtained from both HMM and CRF are provided in Table 5.7. For all four types of observed features reported in the table, the HMM generates higher errors in predicting the labels. The differences in the classification performance of CRF and HMM methods are found to be statistically significant at a 0.05 level. The p‐values of the performed paired t‐tests are less than 0.01 for the different observed features in Table 5.7. It can be noted that adding the tools’ orientation degrades significantly the HMM performance. When using both positions and velocities as observed features, the classification accuracy is improved, but is still lower when compared to the CRF success rates.

The results are consistent with the findings in other studies in the literature (Lafferty et al., 2001; Vail et al., 2007a; Martinez and Kragic, 2008), which report higher classification accuracy of CRF when compared to HMM.