Estimating Motor Control Difficulty in Human–Robot Fine Co-Manipulation Tasks Using Brain Activities

Physical human–robot interaction (pHRI) aims to integrate robotic systems with human cognitive abilities to perform tasks that exceed the capabilities of either entity alone [1]. Two critical aspects must be considered: (1) ensuring safety and reliability through the proper design and control of robotic systems [2,3], and (2) modeling human cognition to predict intentions [4] and responses to robot actions [5,6]. These aspects are deeply interconnected, as physical interactions can be modulated to elicit specific cognitive responses, while cognitive states can reflect the preferred interaction parameters [4,7,8] or modulate feedback from robots to humans [9,10] to enhance perceived safety [11]. For example, in rehabilitation, haptic controllers provide feedback to patients, with their cognitive engagement guiding the adaptation of controller parameters [8,12]. The need for assessing the quality of human–robot interaction using human cognition is even more critical in physical interaction.

When robots physically interact with human operators, it is impossible to model all the human actions to ensure safety and dependability. As a result, compliant control strategies such as admittance and impedance control [13] are utilized to accommodate human operators' natural and unstructured movements [14,15]. These strategies empower a rigid actuator, equipped with either a position or force sensor, to manifest compliant behavior by simulating virtual dynamics via manipulating virtual inertia, damping, or stiffness control parameters. While compliance is crucial to ensuring interaction safety, it alone does not suffice for achieving an intuitive pHRI. Alongside compliance, adaptability to diverse tasks becomes imperative. Studies indicate that controllers adjusting their parameters in accordance with task requirements and human intent outperform constant impedance/admittance controllers. Nonetheless, incorrect control parameter settings can precipitate unstable interactions, thereby compromising operator safety and the system functionality [16].

To minimize effort and maintain stability, the interaction should be continuously assessed, and control parameters should be adjusted in real-time [3,17,18]. Several studies have focused on real-time stability prediction using position [19] or velocity profile analysis [20,21]. Alternative approaches based on frequency-domain analysis have extended stability assessment to interaction forces. For example, Balasubramanian et al. [22] introduced the concept of “spectral arc length,” which evaluates interaction smoothness by analyzing Fourier spectrum variations in speed profiles. While this metric accommodates complex movements, its robustness in co-manipulation scenarios is significantly influenced by individual arm stiffness and task-specific characteristics.

A more reliable approach was proposed by Ryu et al. [23] and later adapted for admittance control by Dimeas and Aspragathos [24]. This method distinguishes between low-frequency components of intended human motion and high-frequency oscillations caused by instability, forming an instability index. We leverage this frequency-based interaction instability index (⁠$IR$⁠) to account for individual variations and user intent.

Instability undermines the efficacy of human-coupled robotic systems and poses risks to user well-being [17,25]. Therefore, investigating the relationship between the cognitive and physical facets of pHRI is crucial for achieving a safer and more organic interaction. The human reaction time to tactile stimulation associated with instability typically ranges from 150 ms to 400 ms [26]. When additional cognitive processes such as motion planning and visual perception are involved (e.g., steering a vehicle), this response time can extend up to 700 ms [27]. Hence, the fast detection time (or even prediction) of any instability in human–robot interaction can help the real-time adjustment of robot behavior to ensure safer physical interaction. In this regard, there have been recent studies that indicate a potential correlation between electroencephalography (EEG) activities and haptic feedbacks in physical human–robot interaction [28–30]. Drawing inspiration from EEG studies in driving assistance systems [31–33], where braking intention has been detected from brain activity up to 380 ms before foot engagement with the brake pedal [34]. We hypothesize that instability in human–robot interaction can be detected from EEG data before they manifest in physical interaction (force data). This proactive detection can potentially enable earlier intervention, reducing reliance on delayed physical responses.

This article investigates the effectiveness of brain activity monitoring in predicting the instability index during physical human–robot interaction and examines how early EEG signals can anticipate instability before they become observable in force data. To achieve this, we conducted a pHRI experiment (Fig. 1) in which participants guided a robot under admittance control along a predefined trajectory. By adjusting the admittance control parameters—some of which were intentionally suboptimal for the task—we induced varying levels of interaction stability. Interaction force data were recorded and analyzed to compute the true instability index. Simultaneously, EEG data were collected and used to predict the instability index via a deep convolutional neural network (CNN) and an elastic net regression (ENR) approaches.

Furthermore, to assess the temporal advantage of EEG-based predictions, we incrementally shifted the prediction window by 0.5 s up to 2 s before the instability was reflected in force data. This allowed us to evaluate how early brain activity can signal impending instability compared to interaction dynamics.

In an experiment approved by the Institutional Review Board, 11 subjects (8 males and 3 females; all right-handed) were recruited from the University of Buffalo School of Engineering, Buffalo, NY. Participants' age ranged from 23 to 34 years (mean = 24, SD = 2.1). The participants were asked to hold the end-effector of a robot with a force–torque sensor and guide it in between a star-shaped trajectory as shown in Fig. 2.

The robot operates under an admittance control framework in which human intent is detected using a force–torque sensor. The robot's end-effector position or velocity is determined based on the magnitude and the direction of the force exerted by the human operator. The robot's response is governed by a virtual dynamic model, which includes parameters such as mass, damping, and stiffness. By adjusting these parameters, the system's sensitivity, stability, and the level of human effort required for interaction can be optimized [16,28,35].

In this equation, $Md$ and $Cd$ represent the virtual mass (inertia) and virtual damping coefficient matrices (set as diagonal matrices), respectively, which are defined for each joint of the robotic system. Here, the diagonal elements are set as $Md$ = [3, 3, 3, 0.1, 0.1, 0.1]. $Cd$ is used to adjust the robot's response by transitioning between high and low damping where an increase enhances motion stability but reduces human comfort, and vice versa. Refer to Ref. [17] for more details on the effect of virtual damping on human effort and motion stability.

The experiment involved altering the controller dynamics of the robot by adjusting the damping coefficient (⁠$Cd$⁠) in the admittance control while keeping other parameters constant. Participants completed separate trials under low-damping (high admittance) and high-damping (low admittance) settings. The high-damping parameters are set as $Cd$ = Diag([80, 80, 80, 10, 10, 10]) and the low-damping parameters are set as $Cd$ = Diag([20, 20, 20, 6, 6, 6]). Low-damping conditions reduced the resistance from the robot, consequently demanding less effort from the human operator to maneuver the end effector. Conversely, high damping increased resistance, elevating the effort required for manipulation. However, under low damping, the robot exhibited heightened sensitivity to force variations, potentially leading to instability during co-manipulation.

To assess the stability of human–robot interactions, we utilized both EEG data and interaction forces between the participant and the robot. The EEG signals were collected using a B-Alert X20 wireless headset (Advanced Brain Monitoring, Carlsberg, CA) at 20 locations (Fp1, Fp2, F7, F3, Fz, F4, F8, T3, C3, Cz, C4, T4, T5, P3, Pz, P4, T6, O1, POz, and O2) and referenced with regard to electrodes located behind each ear on the mastoid bone. The interaction forces were measured using a six-axis force–torque sensor at the robot's end effector.

EEG signals were preprocessed to remove electrical noises by applying band-pass filtering to the EEG signals within the frequency range of 0.3–50 Hz. Subsequently, the EEG signals are epoched into 2-second segments and time-synchronized with force data; any identified bad epochs are excluded from further analysis. We employ independent component analysis with the Picard algorithm to eliminate eye blinks and muscle artifacts utilizing the mne-python library [36].

To quantify the instability of the human–robot interaction, we propose a novel, dimensionless instability index to quantify both low-frequency intended human arm movements and high-frequency unstable oscillations in robots. The premise relies on the distinction in frequency bandwidths between human–arm movements, which typically operate within a relatively low-frequency range [37], and the unstable oscillations in robots resulting from inadequately set admittance levels.

Researchers have extensively explored frequency-domain analysis of end-effector position and interaction force to identify and mitigate unstable oscillations in haptic devices [23,24]. However, the suitability of this analysis depends on the specific type of robot controller being utilized. For instance, Dimeas and Aspragathos [24] demonstrated that spectral analysis of position data remains largely effective for back-drivable haptic devices. In contrast, spectral analysis of interaction forces proves to be more pertinent for robots controlled by admittance. This preference arises because an admittance controller functions like a low-pass filter, attenuating high-frequency force components.

Here $Pf(ω)$ is the fast Fourier transform (FFT) of interaction forces at a frequency of $ω$⁠, $ω0$ is the smallest non-zero frequency component of the FFT, $ωm$ is the maximum frequency component of interest for summation, which must be less than Nyquist frequency, $ωc$ is the cross-over frequency to distinguish between intended and unintended motions which is between $ω0$ and $ωm$⁠. For FFT, the force data were down-sampled to 128 Hz and then segmented to a 2-second data epoch (aligned with EEG signals). A Hamming window with 50% overlap was applied to each segment. Periodograms were averaged to estimate spectral power, and results were normalized according to the overall window power.

The instability index (⁠$IR$⁠) ranges from 0 to 1, with higher values representing higher instability and lower values representing lower interaction instability. The proposed instability index (⁠$IR$⁠) has two parameters namely $ωm$ and $ωc$⁠. The $ωm$ should be chosen so that the effect of sensor noise is minimal and around the maximum range of robot motion. However, the $ωc$ value is intricate, as it is influenced by individual factors such as the user's arm impedance, skill level, and task type. Consequently, diverse users may exhibit distinct $IR$ values for a given task due to variations in $ωc$⁠. To compensate for this individual difference, in this work, the mean value of $IR$ indices is mapped to zero for each subject.

The overall pipeline for predicting the instability index from EEG signals is shown in Fig. 4. We use two different regression approaches: one based on feature extraction (ENR) and the other on raw data (CNN). Each method has strengths and assumptions, influencing how the EEG data are processed. In the first approach, we use ENR, which involves extracting 120 spectral power features from 2-second EEG epochs. We apply a logarithmic transformation to normalize the spectral data, as suggested by Gasser et al. [38]. We also standardize the spectral power values by converting them to z-scores based on each subject's distribution. This approach focuses on deriving meaningful, hand-crafted features from the EEG signals for regression.

In contrast, the second approach uses a CNN that directly processes raw EEG epochs without any feature extraction. The CNN automatically learns relevant patterns from the raw data, avoiding the need for manual feature engineering.

We use leave-one-subject-out cross-validation to evaluate the performance of both approaches. In each fold, the model is trained on all subjects except the tested one (see Fig. 4). For elastic net, a tenfold cross-validation tunes parameters via grid search, while for CNN, a fivefold cross-validation is used due to higher computational demands of deep learning models.

In Eq. (4), the scalar weight $λ$ and the parameter $α$ in Eq. (5) (ranging from 0 to 1) balance the penalty terms, offering a continuum between the L1 norm and squared L2 norm of the vector $w$ [39]. When $α$ tends toward 0, the elastic net method resembles ridge regression, whereas an $α$ value of 1 reflects lasso regression. Solving this regularized optimization problem demands numerical optimization techniques. Ridge regression induces coefficient shrinkage toward zero but does not eliminate any coefficients, retaining all features in the model. Conversely, lasso regression achieves both coefficient shrinkage and feature selection simultaneously. Consequently, some features are assigned zero weights, effectively removing them from the model. The model was trained to predict instability based on the 120 spectral power features derived from the EEG signals.

Alternatively, we also utilized a deep learning approach with a CNN to predict the instability index. CNNs have shown great success in processing time-series data like EEG signals due to their ability to learn spatial and temporal features automatically. The CNN architecture illustrated in Fig. 5 draws inspiration from Refs. [40] and [41].

Next, the data are passed through convolution layers. The first convolution (temporal) uses a 1 × 10 kernel with a stride of 1, producing three channels (using three filters). The second convolution applies a 20 × 20 kernel to capture spatial correlations across electrodes, also with a stride of 1. The output is then squared element-wise, followed by batch normalization and average pooling with a 1 × 50 kernel and a stride of 1 × 10. The pooled output undergoes a logarithmic transformation before passing through another convolution and two fully connected layers with a tanh activation function. The overall architecture is illustrated in Fig. 6. We trained the CNN using the adam optimizer, with an average squared error and a learning rate of $1×10−3$⁠. Training is done for 150 epochs using Nvidia Titan Xp GPU and the PyTorch library.

The instability index $IR$ was predicted using EEG data, and we tested different overlap conditions between EEG epochs to assess whether we could predict the instability before it happens. Specifically, we explored temporal alignments ranging from 0% to 100% overlap between EEG epochs and instability occurrences (Fig. 7).

We tested five overlap conditions, increasing from 0% to 100% in 25% increments. A 0% overlap corresponds to an early prediction of instability by 2 s relative to force data. As the overlap increases (e.g., 25%, 50%, 75%), the prediction becomes less early, gradually approaching real-time detection. At 100% overlap, the focus shifts entirely to correlation analysis rather than early prediction.

The result of frequency analysis of the interaction forces and the individualized instability index (⁠$IR$⁠) are presented in this section followed by the outcomes of the regression analysis for predicting $IR$ from EEG.

The average power spectral density (PSD) of interaction force across subjects is shown in Fig. 8. The power of interaction force is significantly higher in the high-damping compared to the low-damping condition, indicating higher interaction force (human effort). This distinction is visually evident in the vertical shift of the PSD in Fig. 8(a). However, this difference may introduce a potential bias in frequency analysis and instability index calculation. To address this issue, we propose scaling the PSD by the DC gain $(Pf(ω=0))$⁠, as illustrated in Fig. 8(b). This adjustment is essential for mitigating the impact of varying damping conditions on the observed PSD and instability index calculation. Figure 8(b) further reveals that the low-damping condition is associated with higher frequency components beyond human voluntary motion, which corresponds to instability in the interaction force.

To determine the instability index (⁠$IR$⁠), the maximum frequency (⁠$ωm$⁠) is set at 10 Hz, based on the averaged power spectral results. However, individualizing the cross-over frequency (⁠$ωc$⁠) is necessary to account for inter-subject variations in arm impedance and motor control characteristics. This customization is achieved through statistical analysis under both high and low damping conditions. For each subject, the instability index is computed by varying $ωc$ from 0.5 Hz to 4 Hz in 0.25 Hz increments. At each step, the divergence in $IR$ between high and low damping is quantified using F-statistics from analysis of variance (ANOVA). This analysis reveals the trend of F-statistics concerning $ωc$⁠, with the optimal $ωc$ identified at the point where the difference in $IR$ between high and low damping conditions is maximized. This process is shown in Fig. 9. Optimal $ωc$ values, exemplified in Fig. 10 for subjects 5 and 6, are obtained where the F-statistics peak. The optimal cross-over frequencies (⁠$ωc$⁠) of all subjects are listed in Table 1.

This approach requires a significant difference between high- and low-damping conditions to estimate $ωc$⁠. To verify that, we performed a paired t-test on the averaged (across subjects) $IR$ index, which revealed a strong statistical difference (t(11) = −15.1, p < $10−7$⁠) indicating that $IR$ was indeed higher in the low-damping than high-damping condition, validating this index as a quantification of motor control difficulty.

While individualized $ωc$ enhances $IR$ index estimation, it concurrently amplifies between-subject differences. Subjects with lower $ωc$ exhibit larger $IR$ indices, whereas those with higher $ωc$ demonstrate smaller $IR$ indices. Thus, individualized $ωc$ is more appropriate for within-subject analysis rather than between-subject comparisons. To mitigate this issue, an averaged $ωc$ (Table 1) is employed for subsequent analysis. Notably, this average $ωc$ value (2.35 Hz) aligns closely with findings from studies on tremor identification [42,43].

For each subject, the $IR$ index is estimated from EEG data across approximately 600 epochs using CNN and ENR. The distribution of the predicted index for four subjects is shown in Fig. 11.

Despite the differences between actual and estimated $IR$ indices, the CNN model demonstrates reliable performance, particularly in distinguishing between low- and high-damping conditions (p-value < $10−5$⁠). The root mean square error (RMSE) of the estimated $IR$ is averaged for each subject (S1–S11) and compared across different overlap percentages. Table 2 presents this comparison for both approaches. Notably, CNN exhibits superior performance on average compared to elastic net regression. To assess the statistical significance of this improvement, a Wilcoxon signed-rank test is conducted, with results detailed in Table 3.

The results in Table 3 indicate that CNN has notably improved the prediction of $IR$ by up to ∼10.45% compared to ENR. However, this improvement is statistically significant only at 0% and 25% overlap conditions (p-value < 0.05) and not at higher overlap levels. To further assess the magnitude of this improvement, we perform Cohen's d test, which measures the effect size—a standardized measure of the difference between two means. Effect size quantifies the practical significance of the difference, independent of sample size.

In this analysis, 0% overlap has the largest effect size (d = 0.406), indicating a moderate difference between ENR and CNN. The 25% and 50% overlaps show smaller yet noticeable effect sizes (d = 0.33 and d = 0.32, respectively), while 75% and 100% overlap have the smallest effect sizes (d = 0.16 and d = 0.12), suggesting minimal differences between the two models at higher overlap levels. This analysis confirms that the largest differences between the two models occur at lower overlap conditions (0% and 25%), aligning with the Wilcoxon test results, which found significant differences at these levels. This suggests that CNN is a more reliable model for early prediction of instability. However, for real-time detection (higher overlap between force and EEG data), there is no significant difference between the two models.

The stability of co-manipulation depends on the proper adjustment of the virtual damping coefficient in an admittance control strategy. In this study, we explored the feasibility of predicting the instability index from brain activity to enhance adaptive control. As a reference measure, we employed frequency-domain analysis of force data using a 2-second moving window. To predict the instability index, we applied: (1) ENR on power spectral density and connectivity features of EEG data and (2) a deep CNN trained on artifact-free EEG data. Both models demonstrated a reliable estimation performance; however, the estimation quality varied between subjects. The main factor contributing to the varied estimation quality is the functional differences in the subject's motor control skills, which influences their brain activity. This issue can be potentially overcome by training an individualized regression model.

It should also be noted, while the $IR$ index serves as a measure of instability, the precise prediction of its numerical value is less critical in practice than the ability to reliably distinguish between stable and unstable states. Such classification would significantly aid the real-time adjustment of admittance control parameters. Notably, the proposed CNN-based approach successfully and significantly distinguished between stable and unstable conditions, demonstrating its potential for enhancing adaptive human–robot interaction. In future work, we plan to develop an individualized classifier to detect unstable modes from the predicted $IR$⁠, enabling real-time admittance control adjustments. A promising direction for improving robustness is leveraging an ensemble approach that combines ENR and CNN predictions, integrating their strengths to enhance stability detection and control adaptability.

Furthermore, to evaluate the temporal advantage of EEG-based predictions, we shifted the prediction window up to 2 s before instability was reflected in the force data. Notably, the CNN model maintained its performance even with the maximum shift, whereas the elastic net regression model showed a significant decline in performance beyond a 250 ms shift. This finding highlights the potential of EEG-based instability prediction as a proactive measure in physical human–robot interaction. However, it is important to acknowledge the limitations of this study, which involved only 11 subjects performing fine manipulation tasks within a single session. Further research is needed to assess the impact of within-session variability and validate these findings in a larger, more diverse population.

There are no conflicts of interest.

The datasets generated and supporting the findings of this article are obtainable from the corresponding author upon reasonable request.