2. Related work

Recent work has investigated intelligent, adaptive control frameworks that integrate learning-based decision-making for safe human–robot collaboration. For example, Ref. [Reference Jayesh, Altnji, Thiyagarajan, Nanda, Biswas and Chowdhury12] combined Artificial Potential Fields with deep Q-networks for human-aware task planning, enabling dynamic obstacle avoidance and cooperative motion planning, and motivating predictive, context-aware strategies such as MPC in BMI-driven robotics. Similarly, Ref. [Reference Rao and Chowdhury13] developed an audio-aware mobile robot navigation system that exploits environmental sounds for improved decision-making in complex indoor environments. While these approaches rely on external auditory perception, the present work advances auditory integration to the neuro-physiological level by decoding ASSR from EEG for direct robot trajectory control.

SSVEP and ASSR-based BMI for robot control: SSVEP-based BMIs enable fast, high-SNR decoding with intuitive frequency command mapping, and advances in CCA/FBCCA, TRCA, dynamic windows, and LED paradigms [Reference Wong, Wang, Rosa, Chen, Jung, Hu and Wan14] have supported robust online control of robotic arms and mobile platforms. Hybrid SSVEP-EOG systems and SSVEP-driven manipulators demonstrate reliable multi-DOF control. ASSR provides a complementary, eyes-free auditory channel via phase-locked neural responses, reducing visual fatigue. Although standalone ASSR-based robot control remains limited, recent studies validate ASSR paradigms and hybrid audio–visual systems where ASSR serves as a supplementary command source [Reference Wang15].

Hybrid SSVEP and ASSR multi-modal BMI: To mitigate limitations of single-modality BMIs, hybrid systems combining SSVEP with EOG, P300, MI [Reference Tian, Xu, Han, Du, Li, Chen, Xie, Wang, Jiang, Guo, Zhang and Wu16], or auditory steady-state inputs have gained momentum. These approaches improve information transfer rates, robustness, and flexibility, for example, by using SSVEP for mode selection and another modality for trajectory parameters. Recent surveys outline hybrid BMI taxonomies, design challenges, and trends toward adaptive and asynchronous interfaces for assistive robotics [Reference Zhu17]. Although, high SSVEP classification accuracies have been reported in laboratory-based BMI studies, most prior work emphasizes standalone decoding performance rather than integration within safety-constrained robotic control frameworks, often relying on open-loop or simple feedback strategies. Moreover, many BMI-driven robotic systems adopt single-modality or continuous neural control approaches, with limited exploration of discrete shared-control architectures that separate intention decoding from autonomous trajectory generation under model-based predictive control. While ErrP has been investigated for corrective feedback, its integration as a real-time supervisory signal within NMPC-based robotic manipulation remains limited.

Error-related potentials (ErrP) as closed-loop corrective feedback: ErrPs allow online detection of perceived control errors, enabling corrective actions or adaptive relearning [Reference Fidêncio and Klaes18]. When integrated with SSVEP or ASSR decoding, ErrP-based veto or correction improves system reliability, though real-time single-trial classification remains challenging. Ongoing research focuses on generic classifiers, online adaptation, and combining ErrP with reinforcement learning [Reference Salazar-Gomez, DelPreto, Gil, Guenther and Rus19], highlighting its importance for safety and user trust.

Controller implementations for robotic manipulators: Many BMI–robot systems still rely on open-loop mappings or simple feedback controllers [Reference Zeng, Mohammad, Madrigal, Axinte and Keedwell20], whereas model-based methods, particularly MPC, offer predictive control, constraint handling, and optimal input generation. MPC and NMPC have been extensively applied to trajectory tracking and safe human–robot collaboration [Reference Rubagotti, Raimondo, Ferrara and Magni21], with recent work emphasizing real-time feasibility and safety constraints. For instance, [Reference Mirzaeinejad and Shafei22] developed closed-form, prediction-based trajectory tracking controllers for non-holonomic wheeled mobile robots using a recursive Gibbs-Appell formulation, achieving tracking performance comparable to sliding-mode control without requiring online optimization, illustrating the trade-off between computational efficiency and the predictive, constraint-handling capability offered by online NMPC. Studies on UR-series manipulators show MPC consistently outperforming PID in tracking accuracy and constraint satisfaction. Advanced variants, including scenario-based MPC and chance-constrained formulations, provide stronger safety guarantees, which are especially critical for BMI-driven control under decoding uncertainty, delays, and workspace limitations. Recent developments in safety-critical predictive control further strengthen this direction. For instance, [Reference Li, Bi, Li and Gan23] proposed a safety-aware MPC framework for human-in-the-loop intelligent transportation systems, emphasizing constraint tightening and predictive safety margins under uncertain human behavior. Similarly, [Reference Abdelhamid and Mohamed24] developed a human-centered model predictive control architecture incorporating adaptive constraint handling and shared autonomy principles. These works demonstrate how predictive control can formally guarantee safety in collaborative human–machine systems. However, unlike these approaches, which assume continuous human input or behavioral modeling, the present work integrates discrete EEG-decoded task intent with nonlinear MPC while explicitly accounting for neural decoding uncertainty and ErrP-based supervisory correction. This distinguishes the proposed architecture from purely continuous shared-control MPC frameworks. Recent studies have explored advanced control and intelligent optimization frameworks for complex nonlinear systems [Reference Abdelhamid and Mohamed24–Reference Salari, Mirzaeinejad and Mahani28]. These approaches emphasize robustness, adaptability, and optimization under uncertainty. However, they do not consider the integration of neural intention decoding or human-in-the-loop supervision, which is central to BMI-driven robotic systems.

Shared autonomy perspective in EEG–robot interaction: EEG-driven robotic systems can be broadly categorized into (i) continuous neural control architectures, (ii) shared-control frameworks with blended autonomy, and (iii) supervisory discrete-intent paradigms. As discussed in recent human–robot collaboration studies [Reference Li, Bi, Li and Gan23], shared autonomy systems benefit from clearly structured role allocation between human intent generation and autonomous trajectory optimization. The proposed framework adopts this third category, namely a discrete-intent supervisory architecture. Here, the user specifies task-level intent via hybrid SSVEP–ASSR decoding, while the NMPC controller autonomously generates safe trajectories subject to dynamic and workspace constraints. ErrP signals act as a high-level corrective supervisory channel, allowing the human to intervene when perceived execution errors occur. This architecture reduces cognitive demand compared to continuous control while preserving user authority and maintaining formal safety guarantees. As summarized in Table I, prior studies typically address neural decoding (SSVEP, hybrid BMI, or ErrP correction) and predictive control separately. As discussed in the human–robot collaboration studies [Reference Salari, Mirzaeinejad and Mahani28], shared autonomy systems benefit from clearly structured role allocation between human intent generation and autonomous trajectory optimization. Recent safety-critical MPC frameworks in human–machine systems [Reference Li, Bi, Li and Gan23, Reference Li, Li and del Millán29] provide strong theoretical guarantees but do not incorporate EEG-based intention decoding. Conversely, most BMI–robot systems rely on open-loop mappings or low-level feedback controllers without formal constraint-aware predictive control. The proposed framework bridges this gap by integrating hybrid EEG decoding, ErrP-based supervisory feedback, and nonlinear MPC-based constraint-aware trajectory generation within a unified real-time architecture.

3. System architecture

The proposed system architecture integrates EEG-based intention decoding with a constraint-aware robotic execution layer, as illustrated in Figure 2. EEG signals are acquired using either a visual (SSVEP) or auditory (ASSR) stimuli paradigm, which are operated as switchable, mutually exclusive channels for discrete task selection. The proposed framework incorporates two interconnected feedback loops: a control-theoretic loop at the robotic level through NMPC and a human-in-the-loop cognitive feedback mechanism via ErrP signals, enabling supervisory correction of unintended actions.

Table I.

Comparison of representative BMI–robot control frameworks.

Table I Long description

The table compares different BMI-robot control frameworks across several studies. It has six rows and six columns. The columns are labeled Study, Modality, Hybrid, ErrP, MPC, Safety, and Real time. The row labels are Güneşyü & Akin, Zhu et al., Chen et al., Wang et al., MPC-based HRI, and Proposed work. Each row lists the modality used and whether the framework includes hybrid features, ErrP, MPC, safety considerations, and real-time operation. For example, the study by Güneşyü & Akin uses SSVEP modality, does not include hybrid features, ErrP, or MPC, has limited safety, and operates in real-time. The proposed work uses SSVEP + ASSR + ErrP modality, includes hybrid features, ErrP, NMPC, has safety constraints, and operates in real-time.

Figure 2.

Hybrid BMI-NMPC architecture integrating EEG processing, feature extraction, ErrP detection, and robot control.

Following signal pre-processing, feature extraction and classification are performed within the active modality to identify the user-selected command frequency, which is mapped to predefined robotic tasks with offline-generated reference trajectories. ErrP are processed in parallel as a supervisory safety signal; upon detection of execution mismatch or user dissatisfaction, ErrP activation triggers immediate velocity suppression within the control layer. Robot motion execution is governed by a NMPC framework that enforces velocity and safety constraints during trajectory tracking. The NMPC operates independently of the EEG decoding rate, ensuring safe and consistent execution under low-rate and uncertain BMI inputs.

4. Signal pre-processing and feature extraction

The proposed BMI framework uses three modalities of EEG signals, each with a specific role:

• • SSVEP (12, 15, 18) Hz for decoding visual intention,

• • ASSR (35, 40, 42) Hz for auditory-based task selection,

• • ErrP (1–10) Hz for implicit error-related feedback during human–robot collaboration.

EEG signals have very small amplitudes ( $1$ – $100\,\mu$ V), are highly non-stationary, and are easily affected by eye movements, muscle activity, and environmental noise. For this reason, a dedicated pre-processing and feature extraction pipeline is implemented as shown in Algorithm 1 to obtain reliable and robust features for real-time decoding. A detailed step-by-step implementation of the pipeline is provided in the Supplementary Material Section 1.

NOTE: All feature fusion described in this section is performed within a single active modality. The SSVEP and ASSR paradigms are used separately, depending on the task, and no simultaneous audio–visual perceptual fusion is applied.

Let $X_{\mathrm{EEG}} \in \mathbb{R}^{C \times T}$ denote the multichannel EEG signal, where $C$ is the number of EEG channels and $T$ is the number of time samples within a processing window. Similarly, $X_{\mathrm{EOG}} \in \mathbb{R}^{C_e \times T}$ represents the EOG recordings used for artifact removal. The sampling frequency is denoted by $f_s$ . After pre-processing, feature vectors are extracted for steady-state response detection and error-related potential identification. The SSVEP/ASSR feature vector is denoted by $\mathbf{f}_{SS}$ and the ErrP feature vector is denoted by $\mathbf{f}_{ErrP}$ . The final outputs of the pipeline are the estimated stimulus frequency $\hat {f}$ and the ErrP detection decision $d_{\mathrm{ErrP}} \in \{0,1\}$ .

4.1. Pre-processing

4.1.1. Filtering and re-referencing

To remove baseline drift, muscle activity, and line noise, the raw EEG signal $x(t)$ was band-pass filtered between (1–45) Hz and notch-filtered at 50 Hz:

(1) \begin{equation} x_{f}(t) = \left ( H_{bp}(1,45)\cdot H_{\mathrm{notch}}(50) \right ) \ast x(t), \end{equation}

where $x(t)$ is the raw EEG signal, $x_f(t)$ is the filtered signal, $H_{bp}$ and $H_{\mathrm{notch}}$ denote the band-pass and notch filters, respectively, and $\ast$ denotes convolution. All filtering was performed at a sampling rate of 256 Hz. The selected frequency range preserves low-frequency ErrP components as well as frequency-tagged SSVEP and ASSR responses (12 Hz–42 Hz), while suppressing slow drifts and electromyographic noise. The 50 Hz notch filter removes power-line interference common in laboratory environments. The raw EEG signal contains neural activity of interest along with unwanted artifacts such as slow baseline drifts, high-frequency muscle (EMG) noise, and power-line interference. Band-pass filtering between (1–45) Hz preserves the physiologically relevant components for ErrP (low frequencies) as well as SSVEP and ASSR responses (12 Hz–42 Hz), while suppressing slow drifts and high-frequency noise [Reference Diez, Mut, Laciar and Perona10, Reference Wang15, Reference Salazar-Gomez, DelPreto, Gil, Guenther and Rus19]. The 50 Hz notch filter specifically removes line interference common in laboratory environments. This pre-processing ensures that subsequent feature extraction and decoding operate on clean, physiologically meaningful signals, improving both reliability and interpretability of the neural measurements. Representative examples of raw and filtered signals are shown in Figure 3, with modality-specific responses illustrated in Figure 4(a)–(b).

Algorithm 1:

EEG pre-processing and feature extraction pipeline

Algorithm 1: Long description

A textbox containing a step-by-step process for EEG signal processing and feature extraction. The process includes pre-processing steps such as applying band-pass filtering, notch filtering, Common Average Referencing, artifact removal using EOG regression, Independent Component Analysis, and obtaining a cleaned EEG signal. It then describes feature extraction for SSVEP and ASSR, including computing various features and forming a steady-state feature vector. The process also includes estimating stimulus frequency and ErrP feature extraction, forming an ErrP feature vector, and applying adaptive thresholding or classification to obtain a detection decision.

Figure 3.

Representative EEG signals before and after filtering (1–45) Hz, for ErrP (fronto-central), ASSR (temporal), and SSVEP (occipital) modalities. Right: scalp topographies show modality-specific spatial localization.

4.1.2. Common average referencing (CAR)

To reduce global noise and emphasize local neural activity, Common Average Referencing (CAR) was applied after filtering:

(2) \begin{equation} x_i^{\mathrm{CAR}}(t) = x_i(t) - \frac {1}{N}\sum _{j=1}^{N} x_j(t), \end{equation}

where $x_i(t)$ is the filtered EEG signal from channel $i$ , $x_i^{\mathrm{CAR}}(t)$ is the re-referenced signal, $N$ is the total number of EEG channels, and $\frac {1}{N}\sum _{j=1}^{N} x_j(t)$ represents the average signal across all channels at time $t$ . This re-referencing reduces global or common-mode noise shared across electrodes and emphasizes spatially localized neural activity. By subtracting the average signal across all channels, CAR improves the signal-to-noise ratio of the EEG and enhances the stability of subsequent analyses, such as ICA and feature extraction, ensuring that extracted components more accurately reflect cortical activity rather than widespread artifacts. CAR was applied prior to Independent Component Analysis (ICA) to improve the stability of component decomposition and reduce common-mode noise across electrodes.

Figure 4.

(a) SSVEP responses (12,15,18) Hz with the power spectra and PSD; (b) ASSR responses (35,40,42) Hz with spectral comparison.

Figure 4. Long description

Panel A: Three line graphs show SSVEP responses at 12 Hz, 15 Hz, and 18 Hz with amplitude on the vertical axis and time in seconds on the horizontal axis. Below each line graph, a power spectrum graph displays power spectral density (PSD) with frequency in Hertz on the horizontal axis and PSD in microvolts squared per Hertz on the vertical axis. Each power spectrum includes a signal-to-noise ratio (SNR) and phase-locking value (PLV). Panel B: Three line graphs show ASSR responses at 35 Hz, 40 Hz, and 42 Hz with amplitude on the vertical axis and time in seconds on the horizontal axis. Below each line graph, a power spectrum graph displays PSD with frequency in Hertz on the horizontal axis and PSD in microvolts squared per Hertz on the vertical axis. Each power spectrum includes an SNR and PLV. The bottom row features two multi-frequency PSD comparison line graphs. The left graph compares SSVEP frequencies (12 Hz, 15 Hz, 18 Hz) with PSD on the vertical axis and frequency in Hertz on the horizontal axis. The right graph compares ASSR frequencies (35 Hz, 40 Hz, 42 Hz) with PSD on the vertical axis and frequency in Hertz on the horizontal axis.

4.1.3. Artifact reduction

EEG recordings are commonly affected by eye movements, blinks, and muscle activity. To reduce these artifacts, a two-stage cleaning procedure was applied using EOG regression followed by ICA [Reference Zhu17].

EOG regression. Ocular artifacts were attenuated by regressing out EOG contributions from the EEG:

(3) \begin{align} X_{\mathrm{EEG}}^{\mathrm{clean}} = X_{\mathrm{EEG}} - W X_{\mathrm{EOG}}, \qquad W = X_{\mathrm{EEG}} X_{\mathrm{EOG}}^\top (X_{\mathrm{EOG}} X_{\mathrm{EOG}}^\top )^{-1}. \end{align}

where $X_{\mathrm{EEG}} \in \mathbb{R}^{C \times T}$ is the EEG data matrix with $C$ channels and $T$ time samples, $X_{\mathrm{EOG}} \in \mathbb{R}^{C_e \times T}$ is the EOG signal matrix, $W \in \mathbb{R}^{C \times C_e}$ is the regression coefficient matrix, and $X_{\mathrm{EEG}}^{\mathrm{clean}}$ is the EEG signal after removal of ocular artifact components. This step reduces blink and saccade-related activity by removing their linear projection from each EEG channel.

Independent component analysis (ICA). Residual ocular and muscular artifacts were further removed using ICA, which decomposes the EEG into statistically independent components:

(4) \begin{equation} X = AS, \end{equation}

where $X \in \mathbb{R}^{C \times T}$ is the observed EEG data matrix, $A \in \mathbb{R}^{C \times C}$ is the mixing matrix, and $S \in \mathbb{R}^{C \times T}$ contains the independent source components. The independent components were estimated using the unmixing matrix $W_{\mathrm{ICA}}$ :

(5) \begin{equation} S = W_{\mathrm{ICA}} X. \end{equation}

Artifact-related components corresponding to eye blinks, eye movements, and muscle activity were identified based on their temporal and spatial characteristics and removed. The cleaned EEG signal was then reconstructed as

(6) \begin{equation} X_{\mathrm{clean}} = A S_{\mathrm{clean}}. \end{equation}

EOG regression and ICA serve complementary roles in artifact removal. EOG regression specifically targets ocular artifacts by estimating and subtracting the linear contribution of eye blinks and saccades from each EEG channel, thereby reducing their impact while preserving underlying neural signals. ICA further decomposes the EEG into statistically independent components, allowing residual ocular and muscular artifacts to be identified and removed based on their characteristic spatial and temporal patterns. This two-step approach enhances the signal-to-noise ratio of the EEG, ensuring that subsequent feature extraction and decoding rely on physiologically meaningful neural activity rather than contaminating artifacts. Further details on filter implementation, ICA configuration, and artifact rejection criteria are provided in Supplementary Section 1.

4.2. SSVEP and ASSR feature extraction

Let $X \in \mathbb{R}^{N \times T}$ denote the multichannel EEG signal, where $N$ is the number of EEG channels and $T$ is the number of time samples in the analysis window. The sampling frequency is denoted by $F_s$ . The stimulus frequency is denoted by $f$ , and $H$ represents the number of harmonics used in reference signal construction. Filter-bank decomposition produces $K$ sub-band signals $X_k$ . Covariance matrices are denoted by $C$ , and $\bar {C}_f$ represents the class-mean covariance matrix for frequency $f$ . The symbols $\alpha , \beta , \gamma , \delta$ denote weighting exponents used in multi-cue fusion.

4.2.1. Power spectral density (PSD)

Spectral features were computed to capture periodic brain responses associated with SSVEP and ASSR stimulation. The power spectral density (PSD) for each channel $i$ was estimated as:

(7) \begin{equation} P_i(f) = \frac {1}{T F_s} \left | \sum _{n=0}^{T-1} x_i[n] e^{-j 2 \pi f n / F_s} \right |^2, \end{equation}

which quantifies the signal power at frequency $f$ over a finite analysis window where $x_i[n]$ is the EEG signal from channel $i$ , $T$ is the number of samples, $F_s$ is the sampling frequency, $f$ is frequency, and $P_i(f)$ is the power spectral density of channel $i$ . The stimulus-related spectral power was then aggregated across channels and a narrow frequency window $\Delta f$ :

(8) \begin{equation} P(f) = \sum _{i=1}^{N} \sum _{f' \in [f-\Delta f, \, f+\Delta f]} P_i(f'). \end{equation}

where $N$ is the number of EEG channels and $\Delta f$ is the frequency window around the stimulus frequency. To improve robustness against background noise, the normalized signal-to-noise ratio (SNR) was computed as:

(9) \begin{equation} \mathrm{SNR}(f) = \frac {P(f)}{\frac {1}{|\mathcal{N}(f)|} \sum _{f' \in \mathcal{N}(f)} P(f')}, \end{equation}

which emphasizes narrowband stimulus-locked spectral peaks relative to neighboring frequencies. where $\mathcal{N}(f)$ denotes the set of neighboring frequencies used to estimate background noise power. Spectral analysis captures the periodic neural responses elicited by SSVEP and ASSR stimulation. The PSD quantifies the power of each EEG channel at specific frequencies, allowing identification of stimulus-locked responses. Aggregating power across channels and a narrow frequency window emphasizes signals coherent with the stimulation while reducing the influence of unrelated background activity. Normalizing by the power of neighboring frequencies further enhances the signal-to-noise ratio, making the detection of stimulus-specific spectral peaks more robust against spontaneous EEG fluctuations and environmental noise. This approach ensures that the extracted spectral features reliably reflect the neural activity associated with the intended stimuli.

4.2.2. Canonical correlation analysis (CCA)

Canonical Correlation Analysis (CCA) identifies linear combinations of EEG channels that are maximally correlated with frequency-specific reference signals. By exploiting the phase-locked relationship between the brain’s response and the stimulus, CCA effectively enhances stimulus-related neural activity while suppressing unrelated background signals. Including multiple harmonics allows the method to capture nonlinear neural responses that often appear in SSVEP and ASSR paradigms. The stimulus frequency corresponding to the highest canonical correlation is selected as the decoded target, providing a robust and physiologically grounded measure of user intention. It exploits the phase-locked relationship between multi-channel EEG signals $X$ and frequency-specific reference sinusoids $Y_f$ :

(10) \begin{equation} \begin{aligned} Y_f(t) = \begin{bmatrix} {\sin}(2\pi f t) & {\cos}(2\pi f t) & \ \ \cdots \\ {\sin}(2\pi H f t) & \ \ {\cos}(2\pi H f t) \end{bmatrix} \end{aligned} \end{equation}

where harmonics are included to capture nonlinear neural responses, where $H$ is the number of harmonics and $Y_f(t)$ is the reference signal matrix corresponding to stimulus frequency $f$ .

Spatial filters $w_x$ and $w_y$ are obtained by maximizing the canonical correlation:

(11) \begin{equation} \rho (X, Y_f) = \max _{w_x, w_y} \frac {\operatorname {cov}{(w_x^\top X, w_y^\top Y_f)}}{\sqrt {\operatorname {var}{(w_x^\top X)}\operatorname {var}{(w_y^\top Y_f)}}}. \end{equation}

where $w_x$ and $w_y$ are spatial filters and $\rho$ is the canonical correlation coefficient. The frequency yielding the highest correlation indicates the stimulus-locked neural response.

4.2.3. Filter-bank CCA (FBCCA)

To enhance discrimination between closely spaced stimulus frequencies, EEG signals were decomposed into multiple sub-bands $X_k$ [Reference Niu, Zhai and Xia35], and CCA was applied independently to each band:

(12) \begin{equation} \tilde {\rho }_f = \sum _{k=1}^{K} w(k) (\rho _k^f)^2. \end{equation}

where $X_k$ is the EEG signal in the $k$ -th sub-band, $K$ is the number of sub-bands, $\rho _k^f$ is the CCA correlation for sub-band $k$ at frequency $f$ , and $w(k)$ is the sub-band weighting function. This formulation integrates information from fundamental and harmonic frequency components.

4.2.4. Multi-cue fusion

The multi-cue fusion framework integrates complementary information from temporal (canonical correlation), spectral (PSD/SNR), and spatial (Riemannian covariance) domains. The MEC criterion evaluates how well each frequency-specific reference subspace explains the observed EEG, emphasizing stimulus-locked activity. The Riemannian spatial cue captures inter-channel dependencies and provides robustness against noise, electrode shifts, and inter-subject variability by comparing covariance matrices on the manifold of symmetric positive definite matrices. Combining these cues with normalized spectral power via weighted exponents produces a unified decision score, enhancing detection reliability across varying signal-to-noise conditions. This integrative approach ensures that the decoded frequency robustly reflects the user’s neural intent even under realistic EEG variability.

MEC (minimum energy combination). The MEC criterion evaluates the energy of EEG projections onto frequency-specific reference subspaces:

(13) \begin{equation} E_k^{\, f} = \| P_{Y_f} X_k \|_F^2, \quad w_{\mathrm{MEC}}^{\, f} = \frac {\sum _k E_k^{\, f}}{\max _f \sum _k E_k^{\, f}}. \end{equation}

This cue favors frequencies whose modeled responses best explain the observed EEG activity.

Riemannian covariance (spatial cue) [Reference Li, Bi and Yi34]. Spatial dependencies between EEG channels were characterized using covariance matrices on the Riemannian manifold of symmetric positive definite (SPD) matrices [Reference Chevallier, Kalunga, Barthélémy and Monacelli36]:

(14) \begin{equation} \begin{split} C &= \frac {1}{T-1}(X - \bar {X})(X - \bar {X})^\top , \\ d_R(C_1, C_2) &= \left \| \log \!\left (C_1^{-1/2} C_2 C_1^{-1/2}\right ) \right \|_F, \\ s_f &= \exp \!\left (-\gamma d_R^2(C, \bar {C}_f)\right ). \end{split} \end{equation}

where $C$ is the covariance matrix of EEG signals, $\bar {X}$ is the channel mean, $d_R(\cdot ,\cdot )$ is the Riemannian distance between covariance matrices, $\log (\! \cdot \!)$ denotes the matrix logarithm, and $\gamma$ is a scaling parameter. This geometric framework provides robustness to noise, electrode shifts, and inter-subject variability. The Riemannian distance matrix representing inter-condition separability on the SPD manifold is shown in Fig. 1 in the supplementary material.

4.2.5. Final decision rule

Temporal, spectral, and spatial cues were combined into a unified decision score:

(15) \begin{equation} S_f = \frac {(\tilde \rho _f^\alpha )(w_{\mathrm{MEC}}^{\, f}{}^\beta )(s_f^\gamma )(\mathrm{SNR}(f)^\delta )}{\sum _j (\tilde \rho _j^\alpha )(w_{\mathrm{MEC}}^j{}^\beta )(s_j^\gamma )(\mathrm{SNR}(\, j)^\delta )}, \quad \hat f = \arg \max _f S_f. \end{equation}

This fusion balances complementary cues to improve detection robustness across SNR conditions.

Additional implementation details, parameter settings, and supplementary figures illustrating Riemannian distance structures and filter-bank responses are provided in the Supplementary Material Section 1.

4.3. ErrP processing pipeline

Error-related potentials (ErrPs) are low-frequency event-related responses reflecting the neural evaluation of erroneous or unexpected outcomes, most prominently expressed through the N200 and P300 components [Reference Cruz, Pires and Nunes37]. The framework employs a fully analytic and interpretable ErrP detection pipeline designed for real-time operation with minimal computational overhead. Representative ErrP waveforms are shown in Figure 5.

Figure 5.

ErrP components across fronto-central channels highlighting ERN and Pe windows.

4.3.1. Time–frequency features

The time–frequency representation captures transient low-frequency dynamics associated with error-related potentials, while band-specific energy highlights task-relevant spectral modulations. The N200 and P300 ERP components provide robust markers of error processing, and their amplitude contrast serves as a discriminative feature for ErrP detection. Template correlation compares single-trial responses against a subject-specific canonical waveform to quantify similarity with prototypical error responses. By fusing multiple complementary cues, including temporal, spectral, and ERP-based features into a unified ErrP score, the framework enhances detection robustness under realistic EEG variability. Adaptive thresholding allows the system to account for non-stationary signal properties, supporting real-time, calibration-free detection that can respond immediately to user errors during robot interaction. Let $x[n]$ denote the preprocessed EEG signal for a single channel within an analysis window of length $T$ samples, and $F_s$ denote the sampling frequency. Time indices are denoted by $t$ and frequency by $f$ . The symbols $A_{N200}$ and $A_{P300}$ represent ERP component amplitudes computed over predefined temporal windows. The variable $r$ denotes template correlation, and $S_{\mathrm{ErrP}}$ represents the fused ErrP detection score. To capture transient low-frequency dynamics associated with ErrPs, a short-time Fourier transform (STFT) was applied:

(16) \begin{equation} S(t, \, f) = \sum _{n=0}^{T-1} x[n]\, w(t-n) e^{-j 2\pi f n / F_s}. \end{equation}

where $w(\! \cdot \!)$ is a window function, $S(t, \, f)$ is the time–frequency representation, and $T$ is the window length. Band-specific energy was computed as:

(17) \begin{equation} E_{[f_1, \, f_2]}(t) = \sum _{f=f_1}^{f_2} |S(t, \, f)|^2, \end{equation}

capturing task-relevant spectral modulations. where $E_{[f_1, \, f_2]}(t)$ represents the spectral energy within the frequency band $[f_1, \, f_2]$ .

4.3.2. ERP component features

Discriminative ErrP features were extracted from the N200 and P300 components:

(18) \begin{align} A_{N200} &= \frac {1}{t_2 - t_1} \sum _{n=t_1}^{t_2} x[n], \end{align}

(19) \begin{align} A_{P300} &= \frac {1}{t_4 - t_3} \sum _{n=t_3}^{t_4} x[n]. \\[10pt] \nonumber \end{align}

where $[t_1,t_2]$ and $[t_3,t_4]$ denote the temporal windows corresponding to the N200 and P300 components, respectively. The amplitude contrast $(A_{P300}-A_{N200})$ provides a robust marker of error processing.

4.3.3. Template correlation

Single-trial EEG responses were compared against a subject-specific ErrP template using Pearson correlation:

(20) \begin{equation} r = \mathrm{corr}(X_{\mathrm{ErrP}}, \bar X_{\mathrm{ErrP}}). \end{equation}

where $X_{\mathrm{ErrP}}$ is the single-trial EEG segment and $\bar X_{\mathrm{ErrP}}$ is the subject-specific ErrP template.

4.3.4. ErrP score fusion

Multiple complementary cues were combined into a unified ErrP score:

(21) \begin{equation} S_{\mathrm{ErrP}} = \alpha r + \beta (A_{P300}-A_{N200}) + \gamma \frac {E_\theta }{E_\delta + 1}. \end{equation}

where $\alpha , \beta ,$ and $\gamma$ are weighting coefficients, and $E_\theta$ and $E_\delta$ denote spectral energy in the theta and delta bands, respectively.

4.3.5. Adaptive thresholding

To enable online detection under non-stationary conditions, an adaptive threshold was defined as:

(22) \begin{equation} \tau _t = \mu _S(t-L \, : \, t) + \kappa \sigma _S(t-L \, : \, t), \end{equation}

and an ErrP was detected whenever $S_{\mathrm{ErrP}} \gt \tau _t$ where $\mu _S$ and $\sigma _S$ are the mean and standard deviation of the ErrP score over a sliding window of length $L$ , and $\kappa$ is a sensitivity parameter. The choice of a training-free ErrP detection strategy was intentional to minimize calibration effort and improve real-time deployability in human–robot collaboration scenarios. Detailed parameter settings, temporal windows, frequency bands, template construction, and thresholding parameters are provided in the Supplementary Material Section 1 due to space constraints.

5. Human–robot interaction task design

Manipulation tasks were executed using a UR5 robotic manipulator. For each task, reproducible reference trajectories were generated offline using task-space planning combined with inverse kinematics, ensuring consistency across experimental conditions.

The selected tasks correspond to those defined in the $\mathrm{Robothon}^{\circledR }$ 2023 competition task board [Reference So, Sarabakha, Wu, Culha, Abu-Dakka and Haddadin38]:

• • Task 1: Moving a slider along a socket,

• • Task 2: Opening a box lid,

• • Task 3: Picking and placing a probe.

All reference trajectories are generated offline before the experiment and remain fixed during execution.

5.1. Task-space trajectory specification

Table II lists the variables and symbols used throughout this work. These include task-space and joint-space quantities, trajectory interpolation functions, kinematic Jacobians, and the associated operators used in the trajectory optimization and inverse kinematics procedures.

As summarized in Algorithm2 and Figure 6, the end-effector trajectory $x_{\mathrm{des}}(t) = [p_{\mathrm{des}}(t),\, R_{\mathrm{des}}(t)]^\top$ , $t\in [0,T]$ , is generated from waypoints $\{(\tau _i,p_i,q_i)\}_{i=0}^M$ .

Positions are obtained via a $\mathcal{C}^2$ cubic spline, and orientations are interpolated using Slerp. The cubic position spline was constructed with clamped boundary conditions, enforcing zero initial and terminal velocities and accelerations. This guarantees a $\mathcal{C}^2$ continuous trajectory suitable for smooth robot motion. Linear and angular velocities and accelerations $(v_{\mathrm{des}},\omega _{\mathrm{des}},a_{\mathrm{des}})$ follow analytically from these interpolants. Given the joint configuration $q(t)$ , forward kinematics provide $x(t)=f(q(t))$ , and the geometric Jacobian $J(q)$ maps joint to Cartesian velocities. The desired twist is $\xi _{\mathrm{des}}=[v_{\mathrm{des}},\omega _{\mathrm{des}}]^\top$ , and the corresponding joint velocities and accelerations are obtained via the (damped) least-squares pseudo-inverse:

\begin{align*} \dot q_{\mathrm{ref}} = J^\dagger \xi _{\mathrm{des}}, \qquad \ddot q_{\mathrm{ref}} = J^\dagger (\dot {\xi }_{\mathrm{des}} - \dot J\,\dot q_{\mathrm{ref}}), \end{align*}

where damping is applied near singularities, and $z(t)$ may encode secondary objectives such as joint-limit avoidance.

Table II.

Summary of variables used in task-space to joint-space trajectory generation.

Table II Long description

Panel A: Table with symbols and descriptions. Symbols include variables for time, position, orientation, and trajectory parameters. Descriptions provide details on each symbol’s role in task-space to joint-space trajectory generation. Panel B: Continues with additional symbols and descriptions, covering kinematic Jacobians and operators used in trajectory optimization and inverse kinematics.

The resulting reference trajectory is

(23) \begin{equation} \mathcal{T}_{\mathrm{ref}} = \{ q_{\mathrm{ref}}(t),\,\dot {q}_{\mathrm{ref}}(t),\,\ddot {q}_{\mathrm{ref}}(t) \mid t\in [0,T] \}, \end{equation}

computed under UR5 joint constraints.

Once generated, $\mathcal{T}_{\mathrm{ref}}$ can be executed using different control schemes. These reference trajectories serve as fixed motion plans that are executed by the robot controller once triggered by a discrete BMI command. The complete task-space-to-joint-space trajectory generation algorithm and the corresponding trajectory illustration are presented in the Supplementary Material section 2.

Algorithm 2:

Offline task-space to joint-space trajectory generation (overview)

Algorithm 2: Long description

Panel A: The algorithm begins by interpolating end-effector positions using C2 cubic splines and orientations using spherical linear interpolation. Panel B: For each time step, the desired pose and its derivatives are evaluated, and the desired Cartesian twist is computed. Panel C: The task-space motion is mapped to joint space using inverse kinematics with damping, enforcing joint limits and kinematic feasibility. Panel D: The resulting joint positions, velocities, and accelerations are stored and returned as the joint-space reference trajectory.

Figure 6.

Schematic of the trajectory generator module.

6. Control strategy for task execution

6.1. Open-loop control

Each detected SSVEP/ASSR frequency selects a discrete task whose joint-space trajectory $\{\mathbf{q}_{\mathrm{ref},k}\}_{k=0}^{N}$ is precomputed offline via inverse kinematics. The robot executes this plan in open loop at period $T_s$ by applying successive joint setpoints:

\begin{align*} \mathbf{q}_{\mathrm{cmd},k} = \mathbf{q}_{\mathrm{ref},k}, \qquad \mathbf{q}_{\mathrm{ref},k+1} = \mathbf{q}_{\mathrm{ref},k} + \Delta \mathbf{q}_{\mathrm{ref},k}, \end{align*}

where $\Delta \mathbf{q}_{\mathrm{ref},k}=\mathbf{q}_{\mathrm{ref},k+1}-\mathbf{q}_{\mathrm{ref},k}$ .

No feedback correction is applied, so this baseline is limited by decoding latency, noise, and discretizations.

6.2. Nonlinear PID control

To track the precomputed joint-space reference trajectories, a joint-wise nonlinear PID controller was implemented in discrete time with sampling period $T_s$ . The controller outputs desired joint velocity commands, which are sent to the UR5 robot via the RTDE speedJ interface. Let $\mathbf q = [q_1,\dots ,q_6]^\top$ denote the joint configuration, with forward kinematics $\mathbf x = f(\mathbf q)$ and Jacobian $J(\mathbf q)$ . The joint-space tracking error is defined as

(24) \begin{equation} \mathbf e(t) = \mathbf q_{\mathrm{des}}(t) - \mathbf q(t), \end{equation}

which the controller drives toward zero while respecting actuator limits. A per-joint nonlinear PID structure with adaptive proportional scaling, filtered derivative action, integral anti-windup, and velocity feedforward was employed. The resulting discrete-time control law is

(25) \begin{equation} \mathbf u_k = \operatorname {sat}\!\left ( \mathbf K_{p}^{\mathrm{ad}}(\mathbf e_k)\circ \mathbf e_k + K_i\,\mathbf I_k + K_d\,\mathbf D_k + K_f\,\dot {\mathbf q}_{\mathrm{des},k} \right ), \end{equation}

where $\mathbf I_k$ and $\mathbf D_k$ denote the integral and filtered derivative terms, respectively, and $\operatorname {sat}(\! \cdot \!)$ enforces joint-wise velocity limits [Reference Su, Duan and Zheng39]. The adaptive proportional gain improves responsiveness across varying error magnitudes, while feedforward compensation enhances tracking during fast motions. Controller gains were initialized using Ziegler-Nichols step-response tuning and subsequently refined experimentally to achieve critically damped joint responses while respecting actuator limits. Derivative action was low-pass filtered for noise suppression, integral action was bounded to prevent windup, and a velocity feedforward term improved tracking during fast motions. Detailed discrete-time formulations, gain tuning, saturation bounds, and implementation specifics are provided in the Supplementary Material Section 2. The PID controller, tuned as described in Supplementary File S2, serves as a baseline to benchmark the proposed control strategy. While effective in basic tracking tasks, its performance is limited under dynamic or uncertain conditions, highlighting the advantages of the hybrid/BMI-informed controller.

6.3. Proposed MPC controller

Building on predictive-control approaches [Reference Kalaycioglu and de Ruiter40], we propose a NMPC for the 6-DOF UR5 manipulator, incorporating human-in-the-loop safety via ErrP signals. ErrPs act as supervisory cues that adapt the NMPC cost, promoting safe and smooth behavior without explicit stop commands, as shown in Algorithm3. The system architecture is illustrated in Figure 7.

Algorithm 3:

Task-space NMPC with ErrP adaptation

Algorithm 3: Long description

The textbox outlines an algorithm for task-space Nonlinear Model Predictive Control (NMPC) with Error Potential (ErrP) adaptation. It specifies the required inputs such as desired trajectory, joint limits, velocity limits, control limits, prediction horizon, and weight matrices. The algorithm ensures control input and predicted joint states by initializing joint configuration, measuring joint state and end-effector position, computing task-space tracking error, detecting ErrP signal, updating NMPC weights, formulating NMPC optimization, solving the optimization to obtain the optimal control sequence, and applying the first optimal control input. This process is repeated for each timestep.

Figure 7.

Proposed MPC-based safe control architecture for human–robot collaboration.

Figure 7. Long description

A diagram of a control architecture for human-robot interaction using model predictive control. The diagram includes several key components: ErrP Signal, Desired Trajectory, Forward Kinematics, Error Computation, NMPC Objective Function, Constraints, Mapping, and a robotic arm. The ErrP Signal provides feedback to the system. The Desired Trajectory is processed through Forward Kinematics to compute the error. This error is then used in the NMPC Objective Function, which includes a Prediction Model and constraints. The output is mapped to joint space velocity and sent to the robotic arm. Measured data is fed back into the system to complete the loop.

System model. The NMPC state is defined as $x_k = q_k \in \mathbb{R}^6$ , with task-space velocity input $u_k \in \mathbb{R}^3$ . The joint configuration at discrete time $k$ is

(26) \begin{equation} q_k = [q_{1,k},\dots ,q_{6,k}]^\top \in \mathbb{R}^6, \end{equation}

with task-space position

(27) \begin{equation} p_k = f(q_k) \in \mathbb{R}^3, \end{equation}

and tracking error

(28) \begin{equation} e_k = p_{\mathrm{des},k} - p_k. \end{equation}

Modeling assumption and rationale. The NMPC formulation adopts a velocity-level kinematic model of the UR5 manipulator, which is appropriate for task-space control at moderate speeds and widely used in real-time industrial robot control. This choice is motivated by the hierarchical control architecture of such systems, where low-level joint servo controllers operate at high frequency and internally compensate for dynamic effects such as gravity, friction, and inertia [Reference Siciliano and Khatib41, Reference Spong, Hutchinson and Vidyasagar42]. Accordingly, the NMPC layer generates motion references, while the embedded controller ensures dynamic stabilization. Using a kinematic model reflects a complexity-feasibility trade-off. Incorporating full rigid-body dynamics into the NMPC prediction model would significantly increase computational burden and may compromise real-time performance required for human–robot collaboration. Since the proposed framework operates under safety-constrained velocities and smooth motion profiles, dynamic nonlinearities are limited. Neglecting explicit dynamics may introduce prediction errors under aggressive motions. However, these effects are mitigated by (i) embedded joint-level compensation and (ii) the receding-horizon NMPC structure, which continuously corrects model mismatch using measured state feedback [Reference Su, Duan and Zheng39]. This approach preserves stable and accurate trajectory tracking within the considered collaborative operating regime.

Prediction model. The control input is task-space velocity $u_k \in \mathbb{R}^3$ , mapped to joint velocities via the Jacobian pseudo-inverse:

(29) \begin{equation} \dot {q}_k = J^\dagger (q_k) u_k. \end{equation}

Here, $J^\dagger (\! \cdot \!)$ denotes a damped least-squares pseudo-inverse to ensure numerical stability near kinematic singularities. Joint states are propagated using a fourth-order Runge-Kutta discretizations:

(30) \begin{equation} q_{k+1} = q_k + \Phi (q_k,u_k,\Delta t). \end{equation}

ErrPs are modeled as a binary supervisory signal:

(31) \begin{equation} \xi _k = \begin{cases} 1, & \mathrm{if an ErrP is detected at time } k,\\ 0, & \mathrm{otherwise}. \end{cases} \end{equation}

Adaptive NMPC cost. At each timestep, the NMPC solves for $U=\{u_0,\dots ,u_{N-1}\}$ by minimizing

(32) \begin{equation} J = \sum _{k=0}^{N-1} \Big (e_k^\top Q_k e_k + u_k^\top R_k u_k + \epsilon _k^\top D \epsilon _k \Big ) + e_N^\top Q_f e_N , \end{equation}

where $\epsilon _k$ are slack variables for soft constraints. Unlike hard emergency stops, this mechanism yields a gradual reduction in motion through increased control penalties, preserving optimization feasibility and closed-loop continuity. The weights are adapted based on ErrP detection:

(33) \begin{equation} Q_k = Q_0 + \xi _k Q_{\mathrm{ErrP}}, \quad R_k = R_0 + \xi _k R_{\mathrm{ErrP}}. \end{equation}

When $\xi _k=1$ , increased control penalties reduce tracking aggressiveness, biasing the optimizer towards the minimal task-space velocities. Sustained ErrPs naturally lead to $u_k\approx 0$ , achieving a safe stop.

Constraints. Physical and safety limits are enforced:

(34) \begin{equation} q_{\min } \le q_k \le q_{\max }, \quad \dot q_{\min } \le \dot q_k \le \dot q_{\max }, \quad u_{\min } \le u_k \le u_{\max }. \end{equation}

ErrP activation encourages conservative, constraint-consistent behavior in shared human–robot workspaces. ErrP detection acts as a discrete human safety veto rather than a continuous control input. Upon detection at time $k$ , the signal does not overwrite the control command; it adaptively modifies the NMPC cost via weight scheduling, biasing the optimizer toward conservative solutions while preserving feasibility. A qualitative comparison of the proposed approach with alternative control strategies is provided in Table III, highlighting the trade-offs between model fidelity, computational complexity, and real-time applicability.

The comparison in Table III is based on commonly reported characteristics of control approaches in robotic systems. In particular, dynamic NMPC is known for high modeling fidelity but increased computational cost [Reference Chemori, Kouki and Bouani43], while PID controllers offer simplicity but lack predictive and constraint-handling capabilities. Linear MPC methods require model linearization, which may reduce accuracy under nonlinear conditions [Reference Li, Zhang, Wang and Xie44].

6.4. Performance evaluation

The performance of the proposed NMPC controller for the BMI-driven UR5 robot was evaluated against two baseline configurations: open-loop execution and a conventional nonlinear PID controller. The open-loop setup executes planned trajectories without feedback and serves to assess intrinsic planning and execution accuracy in the presence of modeling errors and disturbances. The nonlinear PID controller represents a standard feedback baseline commonly used in industrial manipulators and is implemented solely to benchmark closed-loop tracking performance. The proposed NMPC is evaluated as the primary control strategy, leveraging predictive modeling, explicit constraint handling, and ErrP-based cost adaptation to enable safe and adaptive human–robot interaction.

Table III.

Comparison of control approaches for robotic manipulation.

Table III Long description

A table comparing different control methods for robotic manipulation. The table has five rows and six columns. The columns are labeled Method, Acc., Cost, Constraint handling, Real time, and Remarks. The rows are labeled PID, Linear MPC, Dynamic NMPC, and Kin. NMPC (Proposed). Row 1: PID, Low, Low, No, High, Simple, no prediction. Row 2: Linear MPC, Medium, Medium, Yes, Medium, Needs linearization. Row 3: Dynamic NMPC, High, High, Yes, Low, Accurate, expensive. Row 4: Kin. NMPC (Proposed), Medium–High, Medium, Yes, High, Real-time trade-off.

7. Experimental setup

The proposed NMPC-based human–robot interaction (HRI) framework was validated through simulation and real-world experiments. The system integrates BMI-driven task selection, trajectory planning, and adaptive robot control for safe and reproducible manipulation. Two healthy right-handed adults ( $n = 2$ ) with normal or corrected vision participated after providing informed consent and reported no neurological or musculoskeletal disorders. Three tasks were evaluated (Figure 8) (1) sliding an object along a socket, (2) opening a box lid, and (3) picking and placing a probe. Each task was repeated five times per participant under three control strategies – open-loop, nonlinear PID, and NMPC – to enable quantitative comparison. EEG was recorded using a 16-channel OpenBCI cap (256 Hz, 10–20 system, AFz ground, Cz reference) with impedance below 10 k $\Omega$ . Channels were selected to capture SSVEP (O1, O2, Oz), ASSR (T7, T8), and ErrP (FCz, Cz) activity. Signals were streamed in real time via the Lab Streaming Layer (LSL).

Figure 8.

Illustration of real-time brain–machine interface (BMI)-based robotic manipulation tasks. Each row represents a different task executed through EEG-based control: (a) slider movement, (b) box opening, and (c) probe placement. The sequence highlights the robot’s actions from the initial home position to task completion, demonstrating precise BMI-driven control over the robotic arm and gripper.

Figure 8. Long description

The image contains three rows of photos, each depicting a different task performed by a robotic arm controlled via a brain-machine interface. Each row shows a sequence of actions from the initial home position to task completion. Panel A: The first row shows the robotic arm moving a slider along a socket. The sequence includes the gripper closing, starting to move the slider, and completing the movement. Panel B: The second row illustrates the robotic arm approaching a box, holding the lid, opening the box, and completing the task. Panel C: The third row demonstrates the robotic arm approaching a probe, picking it up, placing it, and completing the task. Each sequence highlights the precise control achieved through the brain-machine interface.

Figure 9.

(a) Human operator with EEG cap using the SSVEP/ASSR-based GUI, task board, and the UR5 robot; (b) Blender-based simulation environment in an HRI setup; (c) Real-time experimental execution and corresponding PyBullet simulation of BMI-controlled robotic tasks, demonstrating consistent task performance between the physical robot (top) and the virtual environment (bottom).

A GUI as shown in fig. 3 in section 3 in supplementary material, allowed switching between SSVEP and ASSR paradigms, displaying the respective stimulus frequencies. A Python script performed real-time filtering, feature extraction, and classification, with detected frequencies triggering predefined UR5 trajectories $\mathcal{T}_{\mathrm{ref}}$ . In SSVEP experiments (Figure 9 (a), participants sat 50 cm from a flickering display and fixated on the target stimulus for (2–5) s. Three circular markers flickered at 12, 15, and 18 Hz, mapped respectively to slider motion, box opening, and probe placement. In ASSR experiments, amplitude-modulated tones blue (35, 40, 42) Hz; carriers: (500, 700, 900) Hz were delivered binaurally at 65 dB SPL in a sound-attenuated room. Although multiple carrier frequencies were used in stimulus generation, ASSR decoding was performed solely based on EEG responses at the modulation frequency, as carrier frequencies are not represented in scalp EEG. The same task mapping was used for both modalities. Non-harmonic frequencies minimized interference and motion artifacts. Each trial consisted of 2–5 s stimulation followed by 2–3 s rest, with randomized task order to avoid learning effects. The UR5 was controlled via TCP/IP (198.160.100.100) for real-time trajectory execution. The proposed NMPC runs with a prediction horizon of $N=10$ and a control timestep $\Delta t = 50$ ms (20 Hz control frequency), solved using CasADi + IPOPT. On an AMD Ryzen 9 7900 12-core processor with 32 GB RAM, the average solve time is 6.4 ms, and the worst-case solve time is 11.8 ms, remaining well below the control period and ensuring real-time feasibility. Solve times remain well below the control period, ensuring real-time feasibility.

9. Limitations

The present study was conducted with only two participants, which limits the statistical power and generalization of the reported findings. The results should therefore be interpreted as a proof-of-concept validation of the proposed integrated EEG robot control framework rather than as evidence of population-level performance. Inter-subject variability is well known in EEG-based brain–machine interface systems, and larger-scale studies are required to assess robustness across diverse users and recording conditions. Future work will include systematic evaluation with a significantly larger participant cohort, quantitative statistical analysis, and user-centered performance metrics to validate the scalability and practical applicability of the proposed framework.

Additionally, the present evaluation was conducted using predefined tasks within a controlled laboratory environment, enabling systematic assessment of the framework’s safety, real-time performance, and BMI-driven adaptation under repeatable conditions. While these experiments establish proof-of-concept feasibility, they do not fully capture the complexity of real-world human–robot collaboration settings. In particular, practical deployment may involve variable task sequences, unstructured human motion, dynamic obstacles, shared workspaces, and changing environmental conditions that were beyond the scope of the current study. To address this limitation, future work will focus on extending the framework to more realistic collaborative scenarios involving continuous human–robot interaction, online task adaptation, dynamic environment changes, and multi-user evaluations. Such studies will enable a more comprehensive assessment of robustness, safety, workload, and operational effectiveness under less structured conditions, thereby providing stronger evidence of the framework’s broader applicability.