Human-robot collaboration (HRC) demands systems capable of dynamically adapting to human actions while maintaining safety and efficiency. Early approaches focused on optimizing task sequencing (Chen et al., 2013; Rahman et al., 2015) by pre-assigning roles to humans or robots, resulting in rigid workflows. While effective in controlled settings, such methods struggle to accommodate real-world variability in human motion and decision-making. Subsequent work adopted leader-follower paradigms, where robots reactively adjust actions based on predefined human workflows (Cheng et al., 2020; Cheng et al., 2021; Ramachandruni et al., 2023; Giacomuzzo et al., 2024). However, empirical studies reveal that human operators prefer retaining task control while also expecting robots to anticipate their needs proactively (Lasota and Shah, 2015). This necessitates real-time monitoring of human actions to enable predictive assistance—a capability absent in existing task-allocation frameworks. Critically, none of these methods actively monitor human actions during execution, limiting their ability to recognize and synchronize with ongoing activities.

Beyond HRC, a rich body of research has investigated human motion analysis, particularly focusing on action recognition (Mao et al., 2023; Yan et al., 2018; Ray et al., 2025) and motion prediction (Mao et al., 2020; Dang et al., 2021; Chen et al., 2023). These methods rely on estimated human joint positions, obtained from vision or inertial sensors, to classify actions or predict motion trajectories. While highly effective on activity recognition benchmarks and gesture-level prediction tasks, the majority of these approaches are designed for offline analysis rather than real-time deployment. A few exceptions demonstrate online operation (Chi et al., 2025; An et al., 2023), but even these focus primarily on recognizing discrete actions in streaming settings. Consequently, they do not provide continuous estimates of human task progression during execution, which is essential for action completion estimation in collaborative scenarios.

Real-time human progress monitoring remains an underexplored topic. To our knowledge, only two approaches that estimate the human progress at the action level have been proposed: Maderna et al. (2019), Maderna et al. (2020) employ Open-end Dynamic Time Warping (OE-DTW) (Sakoe, 1979) to estimate human progression, while Cheng and Tomizuka (2021) propose a Sigma log-normal model for predicting action completion times. The latter reports superior performance over DTW in their evaluations, however, our analysis reveals that OE-DTW, while providing temporal flexibility through nonlinear alignment, suffers from oversensitivity to trajectory shape variations common in real-world human motions. A follow-up work by the same authors (Leu et al., 2022) applies the same method for task planning in a collaborative assembly application.

Dynamic Time Warping (DTW) (Sakoe and Chiba, 1978) is a well-known method for computing similarity between temporally misaligned sequences. Open-end DTW (Sakoe, 1979) relaxes endpoint constraints, enabling partial matches for causal systems. While Maderna et al. (2019) applied OE-DTW for real-time progress estimation, our experiments show its unsuitability to the variability of human actions with extensive real-user studies. We address this limitation through an open-ended variant of Soft-DTW (Cuturi and Blondel, 2017), which replaces DTW’s hard min operator with a differentiable softmin to mitigate local minima. Though open-ended Soft-DTW has shown promise for offline skeleton-based recognition (Manousaki and Argyros, 2023), its application in real-time human progress monitoring remains unexplored.

A deeper limitation persists: standard DTW variants use Euclidean distance, which prioritizes absolute spatial alignment over shape similarity. This proves problematic when human motions preserve geometric structure but vary in speed or amplitude.

Recent works focus on on task-adaptive time warping (Trigeorgis et al., 2016; Matsuo et al., 2023), particularly for aligning machine learning datasets. Trigeorgis et al. (2016) learns complex non-linear representations of multiple time-series based on canonical correlation analysis, while Matsuo et al. (2023) learns a distance metric by training an attention model. However, these methods require large training datasets and full knowledge of the signals, making them incompatible with open-ended scenarios where future data are unknown.

Correlation Optimized Warping (COW) (Nielsen et al., 1998) offers an alternative by maximizing Pearson correlation between signal segments, to match similar segments in fields like chromatography, proteomics, and seismology. However, COW’s rigid windowing sacrifices DTW’s temporal elasticity (Tomasi et al., 2004). Recently, seismological research (Wang et al., 2023) has combined DTW with a windowed correlation-based distance, implementing an offline, one-dimensional approach using a weighted biased cross-correlation distance with windows centered on each sample. While this marks an initial attempt to integrate correlation analysis with DTW, their method fundamentally differs from our requirements for real-time human monitoring.

Our method addresses these gaps by adapting Soft-DTW for real-time open-ended alignment (OS-DTW_EU), overcoming the practical limitations of OE-DTW. Additionally, we introduce the Windowed-Pearson (WP) distance, which computes local Pearson correlations within sliding windows along the trajectory. Unlike COW’s segment-wise approach, WP integrates shape similarity into the DTW framework, enabling both local and global optimal alignment. This combination (OS-DTW_WP) ensures invariance to absolute position shifts and effectively captures local patterns while preserving temporal flexibility–an essential feature for HRC applications, where humans may perform similar motions with varying locations, speeds, and intensities.

For a vector, sequence, or signal a , we denote by a i the element at index i (with indices starting from 0). For a matrix A , A i , j refers to the element in row i and column j with 0-based indexing. Given a vector a , the subvector from index i to j (inclusive) is denoted by a i : j . Submatrices are defined analogously.

Dynamic Time Warping (DTW) (Sakoe and Chiba, 1978) is an algorithm designed to compute the optimal alignment between two time-dependent sequences that may vary in speed, enabling a flexible, nonlinear temporal mapping. Typically, DTW enforces that the warping path starts at the first index and ends at the last index of both sequences. Moreover, each index in one sequence is matched with one or more indices in the other, subject to continuity and monotonicity constraints that preserve the original temporal order. The algorithm outputs both the temporal alignment (commonly referred to as the warping path), and the alignment cost, also known as the DTW cost.

DTW is inherently an asymmetric algorithm, designating one sequence as the reference, b = [ b 0 , … , b n − 1 ] ∈ R n × d , and the other as the query or test sequence, a = [ a 0 , … , a m − 1 ] ∈ R m × d . In this paper, we use the terms sequences, signals, and trajectories interchangeably. Additionally, the DTW algorithms we treat in this paper are generalized to handle multidimensional signals and designed to output the phase τ = [ τ 0 , … , τ m − 1 ] of the query trajectory with respect to the reference trajectory. The phase of a signal, is defined as the normalized progress along a reference trajectory. Each point a i in the query sequence is matched with a point b j i in the reference. The phase τ i of a at i is then defined as: τ i = j i n − 1 ∈ 0,1 . Moreover, while classical DTW uses a pointwise Euclidean (or squared Euclidean) distance, the reported algorithms generalize to an arbitrary distance function δ ( ⋅ , ⋅ ) .

Dynamic Time Warping involves three key steps:

Distance matrix computation: A matrix D is computed to have the distances between all pairs of points of the reference and query sequences.

The detailed algorithm is reported in the Supplementary Material (see Supplementary Algorithm 1), while a visual representation of the DTW alignment is shown in Figure 1.

Classical DTW requires the reference and entire query sequence to compute the optimal alignment. This requirement makes DTW unsuitable for real-time applications when sequences are streamed, as it would need access to future data points to compute the alignment.

Open-end DTW, first employed by Sakoe (1979), is a variant that relaxes the constraint that the last point of the two sequences should match, and computes the alignment which best matches all of the query with a first section of the reference. To do so, after the computation of the cumulative cost matrix R , the index of the reference sequence point, b j * , matching with the last point of the query, is chosen as the one with the minimal cost, namely, j * = arg min j ∈ 0 , n − 1 R m − 1 , j .

Classical DTW, which we will refer to as offline DTW to distinguish it from Open-end DTW, matches the entire sequences, enabling the direct derivation of the phase from the full alignment, with the ability to reference both past and future data from the sequences. In contrast, Open-end DTW can operate on a partially available query trajectory. Although a recursive process could be employed to estimate the alignment of the query trajectory to the truncated reference, this paper focuses on real-time phase estimation, which makes the recursive step unnecessary. Specifically, the phase at time step i ∈ [ 0 , m − 1 ] is estimated causally using only past and current query data, without requiring future information. The phase estimate at step i is given by: τ i = j * n − 1 .

DTW is effective at aligning signals that vary in speed; however, by penalizing both minor and major time shifts equally, it can sometimes produce unrealistic warping paths. This drawback is even more pronounced in Open-end DTW, where no constraint exists to ensure that the final samples of the two signals match.

For offline DTW, this problem has been addressed using path constraints, such as the Sakoe-Chiba Band (Sakoe and Chiba, 1978) and the Itakura Parallelogram (Itakura, 1975), or weighting schemes like Weighted DTW (Jeong et al., 2011), which penalize warping paths deviating from the diagonal. However, these methods are not directly applicable in an online scenario, as the diagonal is unknown.

Soft-DTW (Cuturi and Blondel, 2017) replaces the minimum operation in the forward recursion of the DTW algorithm with a soft-minimum, making the DTW loss differentiable. Specifically, the min operator in the forward recursion of DTW, see Supplementary Algorithm 1 Step 5, is replaced by: min γ a , b , c = − γ ⁡ log e − a / γ + e − b / γ + e − c / γ , γ > 0 , min a , b , c , γ = 0 . to ensure numerical stability when γ > 0 , the soft-minimum is calculated using the log-sum-exp trick: min γ a , b , c = − γ log e − a − μ / γ + e − b − μ / γ + e − c − μ / γ + μ γ ,  if  γ > 0 where μ = max ( a , b , c ) .

Originally designed for time series averaging and clustering, Soft-DTW introduces a smoothing factor that helps mitigate local minima. In particular, the soft-minimum weighs all possible paths, ensuring that slight distortions do not dominate the final alignment. With γ = 0 , the formulation reduces to the standard minimum operation, whereas γ → ∞ results in a cumulative cost equal to the sum of all costs. This formulation allows Soft-DTW to handle temporal variability more effectively. In fact, Janati et al. (2020) show that the Soft-DTW loss is not invariant to time shifts and grows quadratically with respect to the time shift, making it suitable for open-end signal matching.

For completeness, we report a version of the Open-end Soft-DTW algorithm for causal phase estimation in Algorithm 1.

Algorithm 1

Open-end Soft-DTW.

 Inputs:

While the Euclidean distance remains the default metric to measure sample-wise similarity in Dynamic Time Warping, this metric assumes consistent absolute scaling between signals. Though Opend-end Soft-DTW is effective for signals with consistent absolute magnitudes and baseline positions, its reliance on the Euclidean distance makes it sensitive to vertical offsets and amplitude variations, often producing suboptimal warping paths for signals that share geometric structure but differ in execution scale.

Recent approaches like shapeDTW (Zhao and Itti, 2018) address this limitation by converting raw signals into shape descriptors (e.g., piecewise aggregate approximations or discrete wavelet coefficients) prior to alignment.

Inspired by correlation analysis, our approach adapts this shape-sensitive philosophy by introducing a windowed Pearson distance that locally normalizes amplitude differences during alignment. This creates an online-capable method that directly compares trajectory shapes through local correlation analysis, while maintaining DTW’s temporal elasticity. The combination of windowed normalization with open-end alignment proves particularly effective for human-robot collaboration scenarios, where human motion patterns exhibit consistent geometric features but significant trial-to-trial variability in speed and scale.

We formally define the Windowed-Pearson (WP) distance between two signal samples as: δ WP w a i , b j ≔ ∑ k = 0 d − 1 1 − Cov a i − w + 1 : i , k ,   b j − w + 1 : j , k Var a i − w + 1 : i , k Var b j − w + 1 : j , k . (1) to calculate the distance when i < w + 1 or j < w + 1 , we pad the signals with their initial values.

Note that in the one-dimensional case, this distance reduces to the Pearson distance between two segments a i − w + 1 : i and b j − w + 1 : j . Thus it can be rewritten as s δ WP w a i , b j = ∑ k = 0 d − 1 1 − ρ a i − w + 1 : i , k , b j − w + 1 : j , k where ρ ( ⋅ , ⋅ ) denotes the Pearson correlation coefficient between two segments.

By design, the WP distance is invariant to vertical shifts and can effectively capture local correlations. A small window measures similarity between samples based on fine-grained local patterns, while a larger window captures broader, more extended patterns.

Furthermore, we demonstrate that for one-dimensional signals, this distance is equivalent to employing z-normalization as the mapping function to calculate the shape descriptor on a window w , as defined by Zhao and Itti (2018).

Consider two windowed segments a and b of length w . Applying z-normalization, we obtain: a ̃ = a − a ̄ σ a , b ̃ = b − b ̄ σ b , where a ̄ , b ̄ , and σ a , σ b , denote the means and standard deviations of segments a and b , respectively. For z-normalized signals, it holds that: ∑ i = 0 w − 1 a ̃ i 2 = ∑ i = 0 w − 1 b ̃ i 2 = w , and ρ a , b = 1 w ∑ i = 1 w a ̃ i b ̃ i . (10) therefore the squared Euclidean distance between the two normalized segments becomes: ‖ a ̃ − b ̃ ‖ 2 2 = ∑ i = 1 w a ̃ i 2 + b ̃ i 2 − 2 a ̃ i b ̃ i = 2 w 1 − ρ a , b , where the final equality follows from substituting the two previous equalities.

Thus, for z-normalized segments, minimizing the squared Euclidean distance is equivalent to minimizing the Pearson distance (up to the multiplicative constant 2 w ). This makes the WP distance defined in Equation 1 a suitable mapping function for shapeDTW.

OS-DTW_EU and OS-DTW_WP require the tuning of one and two parameters, respectively. Specifically, the smoothing factor γ ≥ 0 and, for OS-DTW_WP, also the window size w ∈ [ 1,2,3 , …   ] .

Assuming access to a training dataset, these parameters can be optimized to minimize the average mean squared error (MSE) between the estimated phase and a ground truth phase. An effective choice for the ground truth is a linear phase evolution, defined as: τ ̄ i = i m − 1 for i ∈ 0 , m − 1 .

Alternatively, the ground truth phase can be computed using offline DTW methods such as Soft-DTW.

In scenarios where the ultimate goal is to minimize a cost function that depends on the phase, the cost function itself can serve as the optimization objective.

While various optimization methods are applicable, we select Bayesian Optimization (Snoek et al., 2012) as our preferred method, as it requires a small number of evaluations of the cost function.

Open-end Soft-DTW, as described in Algorithm 1, is not directly suitable for real-time applications. To address this limitation, modifications are necessary to handle streaming input signals efficiently and to store information in a manner that ensures constant computational complexity at each step, thereby enabling bounded-time computation. Such an adapted algorithm is presented in Algorithm 2.

When a new sample a i arrives, the algorithm first computes d , the i -th row of the distance matrix D (as defined in Algorithm 1), which represents the distances between the new sample and all reference samples. Subsequently, r is computed, corresponding to the i -th row of the accumulated cost matrix R (as defined in Algorithm 1) to update the warping costs.

Algorithm 2

Online Open-end Soft DTW.

 Inputs:

This real-time version only requires storing the last computed rows of matrices D and R , resulting in a constant O ( n ) space and time complexity. This represents a significant improvement over the O ( m ⋅ n ) complexity of the “offline” version described in Algorithm 1. The overall complexity depends also on the computational cost of the distance function δ , which is O ( d ) for the Euclidean distance (where d denotes the number of dimensions of the signals) and O ( d ⋅ w ) for the WP distance (where w represents the window size). Consequently, the per-step time complexity is O ( n ⋅ d ) for OS-DTW_EU and O ( n ⋅ d ⋅ w ) for OS-DTW_WP.

Although each step has constant computational complexity, an optimized implementation is crucial not only for real-time applications but also for offline scenarios where many long sequences need be processed. For brevity, we describe in detail the optimized implementation of Algorithm 1 and then briefly explain how it can be adapted to the real-time version (Algorithm 2).

First, we note that computing the distance matrix requires m ⋅ n evaluations of δ . These operations are mutually independent and thus inherently parallelizable. However, this parallelism does not extend to the subsequent recursion steps, where code efficiency becomes critical. In particular, pre-compilation of this stage can yield significant computational benefits. Given the simplicity of the recursive operations, such optimization is sufficient to ensure real-time feasibility.

When δ is the Euclidean distance, each evaluation has O ( d ) complexity, where d denotes the number of dimensions. The entire distance matrix D can be computed efficiently through vectorized operations across the n and m dimensions. For the WP distance, each evaluation has O ( d ⋅ w ) complexity (where w represents the window size), however, standard scientific computing libraries like numpy and scipy lack built-in support for vectorized computation of this metric.

Calculating the entire D matrix requires computing the correlation between all possible pairs of windows. To achieve this, we construct matrices A ̃ k ∈ R w × m and B ̃ k ∈ R w × n for each dimension k ∈ [ 0 , d − 1 ] , containing all possible windows of the respective signals. Specifically: A ̃ k = a 0 , k a 1 , k a 2 , k a 3 , k ⋯ a m − 1 , k a 0 , k a 0 , k a 1 , k a 2 , k ⋯ a m − 2 , k a 0 , k a 0 , k a 0 , k a 1 , k ⋯ a m − 3 , k ⋮ ⋮ ⋮ ⋱ ⋱ ⋮ a 0 , k a 0 , k a 0 , k ⋯ a m − w + 2 , k a m − w + 1 , k a 0 , k a 0 , k a 0 , k ⋯ a m − w + 1 , k a m − w , k , and similarly for B ̃ k .

Next, we compute the column-wise covariance C k = Cov ( A ̃ k , B ̃ k ) using an efficient method (e.g. numpy.cov). The resulting matrix C k ∈ R ( m + n ) × ( m + n ) satisfies: C i , j k = Cov a i − w : i , k , b j − w : j , k .

The top-left submatrix (of size m × m ) represents the covariance between columns of A ̃ k . The bottom-right submatrix (of size n × n ) represents the covariance between columns of B ̃ k . The top-right and bottom-left submatrices represent the covariance between columns of A ̃ k and columns of B ̃ k .

The matrix D can then be computed efficiently by performing standard operations on these submatrices.

For the online version reported in Algorithm 2, B ̃ k can be computed offline, and for each new sample a i we compute the column-wise covariance between B ̃ k and the vector a i − w : i , k .

Moreover, to ensure numerical stability we add a small value ( 1 0 − 12 ) to the variances at the denominator of the WP distance in Equation 1.

We present two approaches for predicting the completion time of a human action based on the real-time phase estimate provided by OS-DTW. This phase estimate offers valuable insight into the current progress of an action, allowing us to forecast its eventual completion.

The first approach assumes that the user will maintain their current pace, using the observed execution speed to estimate the remaining time. The second approach relies on a nominal execution pace derived from historical demonstrations. Both methods, calibrated with prior data, can enable systems to estimate the future duration of a human action—a capability critical for applications such as assistive robotics and collaborative tasks that require timely intervention.

We assume access to a reference trajectory and a set of training trajectories, each of which is d -dimensional and sampled at a constant time step T s . For convenience, we define t end ( y ) as the function returning the completion time of a trajectory y , and ϕ y ( t ) as the function returning the estimated phase τ for the same trajectory at time t .

The nominal duration t ̄ is defined as the average duration of all training trajectories and the reference trajectory: t ̄ = E y t end y .

Thus, the best a-priori estimate for the completion time, based solely on prior information and the current time, is: t ̂ o t = max t ̄ , t .

We define the nominal estimation method as t ̂ nom t , τ = t + 1 − τ   t ̄ , and the linear estimation method as t ̂ lin t , τ = t τ , where t ∈ R ≥ 0 is the current time and τ ∈ [ 0,1 ] is the estimated phase.

The nominal method assumes execution progresses at the nominal speed, with the remaining time estimated as ( 1 − τ )   t ̄ . The linear method assumes a constant execution speed, scaling the current time t inversely with the estimated completion percentage τ . The linear method is analogous to that employed by Maderna et al. (2019).

The linear estimation method is highly sensitive to phase miscalculations and can be inaccurate when the available trajectory percentage is insufficient to reliably estimate future execution speed. Conversely, the nominal estimation method does not leverage past execution speed as an informative metric for predicting the completion time. To address these limitations, we derive an optimal switching rule based on the estimated phase and elapsed time to transition from the nominal to the linear estimation method.

We begin by calculating the optimal switching time. The mean absolute estimation error (MAE) at time t for the nominal estimation method is defined as MAE nom t = E y t ̂ nom t , ϕ y t − t end y , and for the linear estimation method as MAE lin t = E y t ̂ lin t , ϕ y t − t end y .

To determine the optimal switching time, we calculate the cumulative nominal costs up to each time t and the cumulative linear costs from each time t up to the maximum final time t max = max y t end ( y ) : C nom t = ∑ k ∈ 0 , t MAE nom k , C lin t = ∑ k ∈ t + T s , t max MAE lin k .

The total cost for switching at time t is given by: C total t = C nom t + C lin t .

Thus, the optimal switching time t * is: t * = arg min t ∈ 0 , t max C total t .

A similar procedure can be applied to estimate the optimal switching phase. We define the MAE for a phase τ using the nominal and linear estimation methods respectively as: MAE nom ′ τ = E y t ̂ nom ϕ y − 1 τ , τ − t end y , MAE lin ′ τ = E y t ̂ lin ϕ y − 1 τ , τ − t end y .

The cumulative costs are then: C nom ′ τ = ∫ 0 τ MAE nom ′ z   d z , C lin ′ τ = ∫ τ 1 MAE lin ′ z   d z .

The total cost for switching at phase τ is given by: C total ′ τ = C nom ′ τ + C lin ′ τ .

The optimal switching phase τ * is: τ * = arg min τ ∈ 0,1 C total ′ τ .

We consider a scenario where a human and a robotic manipulator concurrently perform separate tasks. The robot executes a sequence of M robot-task actions R = { a i R } i = 1 M , which are repeated indefinitely. Simultaneously, the human undertakes a sequence of N human actions, denoted as H = { a j H } j = 1 N . The human requires the robot’s assistance to complete a subset of these actions, referred to as joint actions and represented by J = { a l J } l = 1 L , where J ⊆ H .

To formalize the relationship between joint and human actions, we define the operator α ( ⋅ ) to map the index of a joint action to the corresponding human action index, such that a α ( l ) H = a l J . For simplicity, we assume the last human action is a joint action (i.e., a N H = a L J ), and no two consecutive human actions are joint actions (i.e., if a j H ∈ J , then a j + 1 H ∉ J ).

To assist the human, the robot must first complete its current robot-task action a i R before pausing its ongoing task. Once paused, the robot performs a preparatory action (e.g., repositioning or collecting a tool) to prepare for the joint action. After completing the joint action, the robot executes a homing action before either resuming its task or preparing for the next joint action. The sets of preparatory and homing actions are denoted as { a l P } l = 1 L and { a l E } l = 1 L , respectively. A depiction of the task as a hierarchical state machine is provided in Figure 3.

Additionally, we assume access to a set of Q human demonstrations for each non-joint human action a j H ∈ H \ J . These trajectories, denoted as Y j = { y k j } k = 1 Q , consist of the Cartesian positions of the human hand along the x -, y -, and z -axes.

The objective is to minimize the total idle times for both the robot ( Δ total idle R ) and the human ( Δ total idle H ) . This is achieved by optimizing the cost function: C Δ total idle R , Δ total idle H : = Δ total idle R + λ Δ total idle H , (2) where λ > 0 is a weighting coefficient that balances their relative importance.

The collaboration problem is modeled as a finite-horizon episodic POMDP. In this framework, the robot acts as an agent that makes binary decisions between robot-task actions–whether to assist the human or continue its task–while the human is treated as part of the environment. Formally, the POMDP is defined by the tuple ( S , A , T , R , Ω , O ) , where:

• S is the state space;

Note that policy actions a ∈ A should not be confused with the task actions ( a i R , a l P , a l J , a l E ) defined in the previous section.

Each element of the state space S is defined as s = ( a i R , a j H , a l J , Δ start H , y H , Δ idle R , Δ idle H ) , where:

a i R ∈ R is the last robot-task action;

The transition function T ( s , a , s ′ ) : = P ( s ′ ∣ a , s ) is the probability of the state evolving from s to s ′ = ( a i ′ R , a j ′ H , a l ′ J , Δ start ′ H , y ′ H , Δ idle ′ R , Δ idle ′ H ) . The state variables evolve as follows: a i ′ R = a i + 1 mod M R a = 0 a i R a = 1 a l ′ J = a l J a = 0 a l + 1 J a = 1 , with remaining state variables ( Δ start ′ H , y ′ H , Δ idle ′ R , Δ idle ′ H ) updated based on observed interactions. In the next section we describe a model to simulate the evolutions of these quantities.

The reward directly minimizes the cost defined in Equation 2: R s , a , s ′ : = − Δ idle R − λ Δ idle H .

The observation function is defined as O ( s ) : = ( a i R , a j H , Δ start H , τ H ) , where τ H = ϕ y H ( Δ start H ) ∈ [ 0,1 ] represents the estimated completion percentage of a j H , computed from y H using OS-DTW_WP. For each human-only action a j H , we select the first trajectory y 1 j ∈ Y j in the training set as the reference for OS-DTW_WP. Thus, the observation consists of the last robot-task action, current human action, elapsed time, and phase estimate. The definition of the observation space Ω follows accordingly.

Training an RL policy directly on physical hardware is impractical due to time constraints and the constant requirement of human involvement in the task. To address this, we develop a simulated environment that models the collaborative task defined in Section 2.4.1, leveraging human demonstrations to approximate real-world dynamics. This environment enables efficient training of online and on-policy RL algorithms while preserving the POMDP structure formalized in Section 2.4.2.

To model the collaborative task, we assume the duration of each action follows a Gaussian distribution, and estimate them from demonstration data. Specifically, Δ k X ∼ N ( μ X k , σ X k 2 ) , where X ∈ { H , R , P , E } corresponds to human, robot-task, preparatory, and homing actions, respectively.

At the beginning of each episode, we sample from these distributions the durations human actions { Δ j H } j = 1 N , preparatory actions { Δ l P } l = 1 L , and homing actions { Δ l E } l = 1 L . Then, one trajectory y ̃ j is sampled from the set of demonstrations Y j for each non-joint action a j H .

Moreover, to avoid overfitting on the training data, we linearly rescale the time axis of each trajectory y ̃ j to align with each sampled duration Δ j H . As a result, each new trajectory represents either a compressed or stretched version of an actual demonstration. We found this augmentation essential for ensuring robustness and improving the policy’s generalization capabilities.

By employing these quantities, in addition to the state evolution dynamics described in Section 2.4.2, we model the transitions of the POMDP from s = ( a i R , a j H , a l J , Δ start H , y H , Δ idle R , Δ idle H ) to s ′   =   ( a i ′ R , a j ′ H , a l ′ J , Δ start ′ H , y ′ H , Δ idle ′ R , Δ idle ′ H ) as: a j ′ H = β a j H , Δ start H Δ Δ start ′ H = Δ − Δ start H − ∑ k = j j ′ − 1 Δ k H y ′ H = y ̃ j ′ 0 : Δ start ′ H Δ idle ′ R = 0 a = 0 max 0 , ∑ k = j α l − 1 Δ k H − Δ start H − Δ l P a = 1 Δ idle ′ H = max 0 , Δ R + Δ start H − ∑ k = j α l − 1 Δ k H a = 0 max 0 , Δ l P + Δ start H − ∑ k = j α l − 1 Δ k H a = 1 Δ R ∼ N μ R i ′ , σ R i ′ 2 is the duration of the robot-task a i ′ R .

Δ is the duration of the transition: Δ = Δ R a = 0 Δ l P + Δ α l H + Δ l E a = 1 . β is a function that, given the current human action a j H and its elapsed time Δ start H , returns the ongoing human action after a time Δ , namely, β a j H , Δ start H Δ : = arg min a j ′ H j ′ ≥ j | Δ ≤ ∑ k = j j ′ Δ k H .

Moreover we assume that the human always starts from the first action, while the robot is already in operation. As a result, the initial state is non-deterministic. Specifically, we assume the first human action to start when the robot is performing an action a i R .

We derive the initial robot action by sampling from a generalized Bernoulli distribution D R . Namely, each action a i R ∈ R has a probability: P a i R = μ R i ∑ k = 1 M μ R k .

Then, initialize the state as: a l J = a 1 J a i R ∼ D R Δ start H = U 0 , N μ R i , σ R i 2 a j H = β a 1 H , 0 Δ start H y H = y ̃ j 0 : Δ start H Δ idle R = 0 Δ idle H = max 0 , Δ start H − ∑ k = 1 α l Δ k H , where U denotes the uniform distribution.

To solve the POMDP described in Section 2.4.2 within the simulated environment of Section 2.4.3, we select the Proximal Policy Optimization (PPO) (Schulman et al., 2017) algorithm. Our choice is motivated by four key factors:

Hybrid State Space Handling: PPO natively supports state spaces combining both discrete ( a i R , a j H ) and continuous observations ( Δ start H , τ H ) , eliminating the need for algorithm modifications.

While alternative algorithms could be considered, PPO avoids several pitfalls that make them less suitable for our application. For example, Deep Q-Networks (DQN) (Huang, 2020) are limited to discrete action spaces and would require explicit discretization of continuous inputs, while also struggling with partial observability. Trust Region Policy Optimization (TRPO) (Schulman et al., 2015) shares many of PPO’s theoretical guarantees, yet, it incurs significantly higher computational overhead without providing empirical gains for binary action spaces. In the case of Soft Actor-Critic (SAC) (Masadeh et al., 2019), its design for continuous control introduces unnecessary complexity when applied to discrete action spaces. Overall, PPO emerges as the most suited choice for the PACE framework, striking an optimal balance between simplicity, sample efficiency, and performance.