Multi-query Robotic Manipulator Task Sequencing with
Gromov-Hausdorff Approximations

The RTSP is a well-studied problem, and approaches can generally be divided into offline and online approaches. Offline approaches tend to formulate the problem as a generalised TSP (GTSP) [10] (also known as the set TSP) where task inverse kinematics (IK) solutions are formed into subsets and the shortest tour between subsets must be found [1]. Determining an optimal sequence generally takes several minutes to hours for 6-DOF and higher robot arms, even when ignoring the cost of the low-level motion of the robot and are thus unsuitable for online planning [11, 12, 13].

Online methods sacrifice sequence optimality for faster planning. For example, if the kinematic redundancy of the arm is ignored and only task-space distances are considered, then the sequencing problem can be treated as a TSP [14, 15]. While this results in much faster task sequencing, task-space costs can heavily underestimate the true cost of moving between tasks due to the non-linear nature of high-dimensional robot arms [16].

Decoupled approaches overcome this limitation by first computing a sequence in task space and then determining an optimal goal configuration assignment. For example, RoboTSP [2] initially computes a sequence using task-space distance and then assigns optimal joint configurations to each task using a graph search algorithm. The work in [3] adds an extra filtering step to RoboTSP which removes dissimilar candidate IK solutions from the graph search which results in faster task sequencing yet similar execution times. While these decoupled approaches are an improvement over task-space sequencing alone, the choice of configuration assignment is dependent on the initial task sequence and may still produce highly suboptimal trajectories.

More recently, Cluster-RTSP [4] was proposed to overcome the shortcomings of decoupled approaches by directly sequencing in configuration space. To make the problem tractable, unique configurations are first assigned to each task based on a best-fit heuristic. These configurations are then clustered into groups and a sequence is found by solving inter-cluster and intra-cluster TSPs individually. While this method achieved a large reduction in task execution and computation times compared to RoboTSP, the travel costs used for sequencing were still Euclidean distance between assigned task configurations which may not effectively capture costs incurred by the true low-level motion, for example, due to obstacle avoidance. Furthermore, outlier tasks with highly dissimilar IK solutions can heavily impact the overall quality of the task sequence.

Learning-based approaches have been proposed to solve RTSPs. For example, in [17] the problem was framed as a Markov decision process and a reinforcement learning algorithm learned a policy that outputs a time and energy efficient task sequence and corresponding goal IK solution for a robot spot welding task. While results show an overall reduction in computation time and trajectory efficiency compared to offline methods, they assume no obstacles between tasks. Given the black-box nature of learning-based approaches it is unclear how well they generalise to the more complex cluttered environments we consider in this work.

The Mobile Manipulator RTSP extends this problem to mobile manipulators. The work in [18] approaches this problem by clustering tasks based on the arm’s reachability from a finite set of base poses. A minimum set of base poses is then chosen and RoboTSP is used to sequence each cluster. We show that $\epsilon$-GHAs can be naturally extended to mobile bases and enhance subspace allocation, allowing for greater flexibility in subspace assignments in addition to operating in potentially larger workspaces.

While RTSP methods exist that explicitly consider trajectory costs, they either simplify the problem or are intractable for anything other than a small number of tasks. For example, in [19] configurations are arbitrarily chosen for each task and motion plans and generated between them. In [5] the RTSP was formulated as a combined set covering problem and TSP (SCTSP) and trajectories were exhaustively computed between all possible goal configuration assignments which took 3 hours for 400 tasks. In [20] a genetic algorithm was used to find an optimal sequence that considered obstacle avoidance in its optimisation objective. However, problems with only 15 tasks were considered and computation times were not reported. The work in [21] formulated a robotic laser welding problem as a TSP with neighborhoods (TSPN) and was given a computation budget of 600 seconds for up to 71 tasks; however, obstacles were ignored. In contrast, HAP computes a task-space decomposition that considers robot kinematics and prior obstacle knowledge in a reasonable amount of time and facilitates rapid online multi-query task sequencing with high-quality trajectories.

Our framework has useful properties for motion planning applications. As we will show, trajectories contained within $\epsilon$-GHAs are useful as seed trajectories for trajectory optimisation methods [22]. These methods are well-known to be sensitive to the initial seed trajectory and susceptible to failing due to local minima [23]. Trajectories contained within $\epsilon$-GHAs mitigate this issue since they consider manipulator kinematics and a priori knowledge of the environment.

Work in integrated task and motion planning [24] typically utilise some decomposition of the workspace for problem reduction. For example, in [25] authors characterise a lower-dimensional subspace of valid grasp and place configurations. A roadmap is then constructed over these subspaces to search for transit/ transfer modes. However, the construction of these subspaces is mainly concerned with feasibility rather than efficiency.

Additionally, other aspects of our approach share similarities with trajectory library methods [26, 27, 28], workspace decomposition [29, 30] and manipulator coverage planning [31, 32]. However, our approach differs in that the task-space decomposition and paths through them are constructed specifically in a way to facilitate efficient robot manipulator task sequencing. Furthermore, we provide practical guarantees on the lengths of paths through the decomposition which leads to predictable and well-defined behaviour.

$\mathcal{C}$ represents the configuration space of the robot. In this work we target serial manipulators. The workspace, $W$, is the 3D Euclidean workspace, $W=\mathbb{R}^{3}$. Given a configuration $q\in\mathcal{C}$, $A(q)\subset W$ denotes the space occupied by the robot model at configuration $q$. $\hat{m}$ $\subset\mathbb{R}^{3}$ is an approximate model of the environmental obstacles. We assume access to a collision checking process that reports whether the arm is in collision with $\hat{m}$ or with itself. The obstacle region is defined as $\mathcal{C}_{\text{obs}}=\{q\in\mathcal{C}\mid A(q)\cap$$\hat{m}$ $\neq\emptyset\}$, from which we obtain the free space region $\mathcal{C}_{\text{free}}=\mathcal{C}\setminus\mathcal{C}_{\text{obs}}$. A task, modelled as a 6D pose $t$, typically has a set of inverse kinematics (IK) solutions $\Xi(t)$. Then, the task space $\mathcal{T}\subset SE(3)$ is the set of poses of the robot’s end effector for which valid IK solutions, $Q(t)\subset$ $\Xi(t)$ exist. Any IK solution $q$ is considered valid if $q\in\mathcal{C}_{\text{free}}$. $\hat{T}$ is a discrete approximation of the subset of $\mathcal{T}$ where the robot is expected to operate frequently. Summary of notation can be found in Table LABEL:table:notation.

To complete a given task $t\in\mathcal{T}$ chosen from the task space $SE(3)$ the manipulator must position its end effector at pose $t$ while avoiding collision. We are interested in finding a minimum cost path in $\mathcal{C}_{\text{free}}$ that completes task $t$. The manipulator’s path is modelled as a discrete sequence of configurations $\pi=\{\pi[1],\ldots,\pi[N]\mid\pi[i]\in\mathcal{C}_{\text{free}}\}$. We measure the length of a path using a metric on the configuration space,

It is convenient to consider the starting pose of the end effector to be the goal pose of the previous task. A manipulator may be able to achieve the goal pose with multiple, possibly infinite, configurations. We are therefore interested in the set of paths (sequence of configurations) of varying length $N$ between task $t_{j}$ and $t_{l}$, $\Pi(t_{j},t_{l})=\{\pi\mid\pi[1]\in Q(t_{j}),\pi[N]\in Q(t_{l})\}$, leading to the following problem definition.

The operational scenarios that motivate this work often involve more than one task. The manipulator is given an unordered set of $M$ tasks $T=\{t_{i}\}_{i=1}^{M}\subset\mathcal{T}$. This set of tasks can be viewed as a batch-query scenario where all tasks must be completed while minimising total cost. Thus, it is necessary to choose a sequence of tasks, imposing a total ordering over $T$, in addition to repeatedly solving Problem 1. The robot task sequencing problem can thus be formulated as follows.

Problem 2 does not admit a tractable solution. For each task pose $t_{j}\in T$ to be visited, there can be multiple valid configurations in set $Q(t_{j})$¹¹1For infinite sets, as is the case for redundant arms, we produce a finite set by discretising around one of the free joints. A direct solution requires simultaneously choosing the optimal configuration from $Q(t_{j})$ and sequencing the tasks. In other words, this case contains an instance of the generalised TSP (GTSP) [10] (Fig. 3b), which is even more complex than the standard TSP [33] (Fig. 3a), and does not allow for real-time planning. The number of possible sequences to evaluate is $\mathcal{O}(|Q(t_{j})|^{M}\times M!)$, where $|Q(t_{j})|$ is the cardinality of the largest set of IK solutions for any task $t_{j}\in T$. Furthermore, evaluating the path cost for each configuration pair requires solving a motion planning problem, which is itself PSPACE-complete [34]. These computational complexities, which we wish to avoid, motivate our approach to solving Problem 2.

Reducing the set of IK solutions of each task $t_{j}$ to a unique solution $q_{t_{j}}\in Q(t_{j})$ transforms the GTSP in Problem 2 to a classical TSP, which is easier to solve. We therefore search for a mapping $\theta$ from each task $t_{j}$ in a discrete approximation $\hat{T}$ of the full task space $\mathcal{T}$ to a suitable unique IK solution $\theta(t_{j})\in Q(t_{j})$.

To ensure that task sequencing using this mapping guarantees short, smooth and collision-free paths between tasks in both task and configuration space, the map $\theta:\hat{T}\to\mathcal{C}$ is chosen to be an approximate isometry. Intuitively, an approximate isometry enforces that two tasks close in task space remain close in configuration space after mapping. Once $\theta$ is found, the GTSP in Problem 2 is successfully reduced to the classical TSP problem as desired, and task sequencing can be solved efficiently. As a result, the solution of the reduced problem is a near-optimal solution to Problem 1, producing smooth, short paths between tasks.

However, the kinematics of a robotic manipulator working in the task space may not necessarily admit a single approximate isometry $\theta$. For example, in Fig. 2a there is no single $\theta$ that allows the arm to travel the short task-space distance required without taking a long path through configuration space. To handle such cases, we relax the constraint of enforcing a unique $\theta$ and consider a finite set of maps. Thus, we propose finding $K$ mappings $\{\theta\}_{t=1}^{K}$ from $K$ subspaces in task space to distinct subspaces in configuration space, i.e., $\theta^{i}:\hat{T}^{i}\rightarrow\mathcal{C}^{i}$. Within each subspace, the corresponding approximate isometry $\theta^{i}$ guarantees short, smooth and collision-free paths in both task and configuration space.

This effectively decomposes the task space into subspaces that may overlap, leading to cases where a single task is assigned multiple IK solutions rather than our single desired one. As such, the GTSP problem is no longer reduced immediately to a classical TSP one. To solve Problem 2 we must first disambiguate which IK solution to select in overlapping subspace regions. Then, we may solve a classical TSP problem within each subspace independently. The final reduced problem and approach is illustrated in Fig. 3c.

The offline pre-computation stage begins with computation of the $\epsilon$-GHA(s). We are given an anticipated environment $\hat{m}$ and a set of tasks $\hat{T}$ representative of online scenarios. We assume the environment is not entirely known in advance, but a general model $\hat{m}$ is available that approximates what is expected. For example, $\hat{m}$ might represent a general bookshelf structure including the shelves and case. However, $\hat{m}$ need not include all objects within the bookshelf, as these details may be unknown a priori and discovered later.

An undirected graph $G$ is created with nodes corresponding to poses $t\in\hat{T}$ and edges $\hat{E}$ formed by connecting nodes within a specified radius of each other. An example $\hat{T}$ and $G$ is shown in Figs. 5a-LABEL:sub@fig:offline_edges. An $\epsilon$-GHA is then generated from $G$ according to Alg. 1 and is represented as a subgraph of $G$; that is, $(\theta,E)$ where $E\subseteq\hat{E}$. An example subgraph representation is shown in Fig. 5c.

Algorithm 1 is initialised by assigning all nodes to $\hat{T}_{\text{open}}$. A node stays in $\hat{T}_{\text{open}}$ until a unique IK solution is assigned. While $\hat{T}_{\text{open}}$ is not empty, $\theta$ is found via a generate map algorithm procedure (Alg. 2). Algorithm 2 searches for a $\theta$ that minimizes the objective cost, the sum of all minimum cost paths $\pi^{*}(t_{0},-)$ from a root node $t_{0}\in G$ to all other nodes. That is,

where the cost of any path $\pi(t_{0},t_{N-1})=\{\theta(t_{0}),\ldots,\theta(t_{N-1})\}$ of $N$ nodes is defined as

To ensure that paths have bounded lengths and do not involve large, unnecessary arm movements, $\theta$ is additionally constrained such that it is an approximate isometry, or an $\epsilon$-Gromov-Hausdorff approximation ($\epsilon$-GHA) [35], defined below.

It should be noted that $d_{T}$ and $d_{C}$ operate on different metric spaces. However, this can be accounted for through appropriate choice of epsilon.

Depending on the value of $\epsilon$, topology of $\hat{m}$ and robot kinematic structure, some nodes may still have undefined mappings after a single iteration of the algorithm. It may then be necessary to search for multiple $\epsilon$-GHAs, $\{\theta^{i}\}_{i=1}^{K}$, and corresponding edges, $\{E^{i}\}_{i=1}^{K}$, that map a covering set of subspaces $\{\hat{T}^{i}\}_{i=1}^{K}\subset\hat{T}$ to a set of disjoint configuration subspaces $\{\mathcal{C}^{i}\}_{i=1}^{K}\subset\mathcal{C}$.

Note that poses in the example in Fig. 5a have the same orientation. For task spaces with differing orientations the algorithm runs identically. As long as a valid task-space distance metric is used, (6) holds.

The generate map algorithm as outlined in Alg. 2 is based on Dijkstra’s algorithm and attempts to find a unique IK solution for each task in $\hat{T}$ such that the objective cost in (4) is minimised. The process is visualised in Fig. 6. It begins by assigning an undefined mapping to all nodes except for the root node which is mapped arbitrarily to one of its IK solutions. The rest of the procedure is carried out as in the original Dijkstra’s algorithm with a priority queue $\mathbb{Q}$ [36], with modifications to the node expansion step where we compute the unique IK solution mapping. Here, the set of neighbouring nodes of $t$, $\mathcal{N}_{t}$, is run through the functions shown in Algs. 3 and 4.

In get mapping (Alg. 3), if a node $u\in\mathcal{N}_{t}$ has already been assigned an IK solution $q_{u}$, this solution and the edge cost,

are returned and the $\epsilon$-GHA $\theta$ is not updated. Otherwise, the IK solution that gives the minimum $l(u,t)$ while satisfying the $\epsilon$-GHA constraint in (6) is returned. The returned $q_{u}$ and $l(u,t)$ are passed as a candidate to update (Alg. 4). In update, if the candidate path cost from root node to $u\in\mathcal{N}_{t}$ is less than the current cost then $u$ is mapped to $q_{u}$, the path cost is updated, and an edge is created between $u$ and $t$. Additionally, if $u$ is not in $\mathbb{Q}$, it is added.

The above process is repeated for all IK solutions of the root node. If required, the algorithm is run $K$ times to find multiple $\epsilon$-GHAs. The resulting one or more $\epsilon$-GHA(s) $\theta^{i}$ that minimise $J(\theta^{i},t_{0})$ are returned. All nodes are assigned an IK solution before being placed in the queue. It is also important to note that these assignments do not change during an iteration of the routine, as this could result in unstable behaviour and violation of the $\epsilon$-GHA condition due to changing edge costs. However, path costs may change after finding shorter paths through $G$.

In contrast to the original Dijkstra’s algorithm, all nodes are initialised with a non-infinite path cost $c_{\text{max}}$. Thus, finding a small number of $\epsilon$-GHAs with greater coverage will be favoured, particularly in earlier iterations.

With each offline task mapped to a unique or reduced set of IK solutions via $\theta$ or $\theta^{i}$ found offline, the online planner can solve a classical TSP as a proxy for solving the task sequencing problem in Problem 2. However, the set of online tasks $T$ provided may not necessarily align with the offline set $\hat{T}$ used to generate the subspace mappings $\theta^{i}$ in the offline stage, as illustrated in Fig. 5d. As such, each task $t\in T$ must first be assigned to one of the subspace mappings $\theta^{i}$. This matching process and subsequent task sequencing is outlined in Alg. 5 and visualised in Fig. 7.

Experiments consist of sets of randomised trials where tasks are sampled uniformly at random from the environment’s task space. The number of tasks range from 5 to 30 and for each setting we run 50 trials to gauge the robustness of each method.

We compare HAP with competitive contemporary RTSP solver baselines. Both methods compute a sequence of goal IK solutions and then use the same trajectory adaption process as HAP, except with a straight-line trajectory prior for the trajectory optimisation planner:

• •

  RoboTSP [2]: Initially solves a TSP in task space. Given this sequence, an optimal joint configuration is assigned to each task using a graph search algorithm such that the total path length (Euclidean distance in configuration space) through each configuration is minimised.

• •

  Cluster-RTSP [4]: Begins by assigning a unique configuration to each task based on a best-fit heuristic. This essentially tries to choose a set of configurations that are all close to each other. The configurations are then clustered into similar groups, based on proximity in configuration space. A configuration sequence is then found by solving inter-cluster and intra-cluster TSPs individually.

The following metrics are evaluated:

• •

  Planning success rate refers to the percentage of tasks for which the trajectory adaption process succeeded in motion planning a trajectory in $\mathcal{C}_{\text{free}}$.

• •

  Motion planning time refers to the total computation time taken for the trajectory adaption process to compute motion plans for all tasks.

• •

  Average execution time refers to the average time taken per task for a successful trajectory execution, assuming the arms are operating at their maximum joint velocities. Note that for HAP we count inter-subspace trajectories, that is the transition to and from the home configuration, as one task.

• •

  Maximum trajectory jerk refers to the maximum joint jerk norm of an executed trajectory. This is computed numerically using the finite difference method.

• •

  Sequencing time for HAP refers to the time taken to match tasks, retrieve subspace paths and then compute a sequence (Alg. 5 lines 2-7). Additionally, online task IK solution computation time is included for all methods.

The robot models used are a 7-DOF Rethink Robotics Sawyer and a 6-DOF Universal Robots UR5. The UR5 also has the option to be mobile, which we refer to as HAP-Mobile. We additionally include a variant of HAP which, similar to RoboTSP and Cluster-RTSP, utilises a straight-line trajectory prior for trajectory optimisation, referred to as HAP-no-prior.

All benchmarks are run on the bookshelf environment shown in Fig. 5. The task-space decomposition is carried out with the empty bookshelf and in the online scenarios objects are added. The computed subspaces are shown in Figs. 14-16.

In implementing Alg. 1, the number of root nodes, $\|\hat{T}_{\text{root}}\|=10$. The maximum number of $\epsilon$-GHAs was set to $5$ for the UR5 experiments. Sawyer experiments only utilised a single $\epsilon$-GHA due the flexibility in IK solutions afforded by its kinematic redundancy. In implementing Algs. 2-4, all nodes are initialised with path cost $c_{\text{max}}=5.0$. The subspace exploration penalty $\rho=2.0$. The subspace distance biasing penalty $\rho_{s}=0.02$. Experiments with UR5 model use $\epsilon=0.35$, and those with Sawyer model use $\epsilon=0.85$. Distance $d_{C}$ is defined as $L_{\infty}$, and $d_{T}$ as $L_{2}$. Both metrics use uniformly weighted distances.

For the online planner (Alg. 5), the number of closest neighbours used when retrieving $\pi^{*}(t_{j},t_{l})$ is $k=10$. The threshold used for terminating the search over mappings when matching a task configuration is $L_{2}$ distance $0.7$. Google’s or-tools [38] package is used as the TSP-solver in (3).

The discretised task space, $\hat{T}$, used is shown in Fig. 5a. The home pose used is shown in Fig. 5d. For the mobile base experiments, a discrete set of possible base positions is defined uniformly along the $y$-axis, parallel to the bookshelf, see Fig. 16.

The trajectory optimisation algorithm we use for online adaption is TrajOpt with default implementation from [39]. The fallback probabilistic method used is BIT* with a timeout limit of 2 seconds and otherwise default implementation from [40].

The simulation environment used is OpenRAVE with IKFast kinematics solver [41]. To generate a finite set of IK solutions for the 7-DOF Sawyer, IKFast sets the second DOF as a free joint with discretisation increments of $0.01$ radians.

For the benchmark methods RobotTSP and Cluster-RTSP, the default parameters are used in the provided code implementations [42, 43]. The only difference being that the maximum number of clusters for Cluster-RTSP was modified to be the number of tasks.

Time benchmarks for the Sawyer and UR5 bookshelf experiments are shown in Figs. 8 and 9. Average maximum trajectory jerk values are shown in Tables II and II.

The box-and-whisker plots in Fig. 10 show that HAP variants achieved favourable success rate relative to baseline methods. Notably, the UR5 benefits from subspace trajectory priors more so than the Sawyer which is able to utilise its kinematic redundancy to better navigate around obstacles.

Motion planning times for HAP variants in the UR5 experiments greatly outperform benchmarks with approximately up to 3x speedup compared to robotsp and 2x compared to cluster-RTSP with consistently higher planning success rates and whilst retaining similar task sequencing time scaling. For a plan to fail, both TrajOpt and BIT* must timeout. Thus, low planning success rates explain the higher motion planning times exhibited by the baseline methods.

HAP variants and Cluster-RTSP saw a downward trend in execution times as the number of tasks increased for both experiments, whereas RoboTSP’s fluctuated for the UR5 case. Whilst HAP’s execution times were similar to Cluster-RTSP, and RoboTSP in the case of Sawyer experiments, the average maximum trajectory jerk values in Tables II and II were up to approximately 5x lower with low variance. It should be noted that trajectory execution times describe successful trials only and should be interpreted in conjunction with corresponding planning success rates.

Furthermore, the HAP-mobile variant achieved comparable results to the static case, confirming HAP’s ability to generalise to mobile bases. Interestingly, the mobile and static variants utilised on average 2.45 and 1.9 subspaces per trial, respectively. This can be explained by the greater subspace coverage and diversity provided by the mobile base (see Fig. 15), hence the slightly better planning and success rates.

Total offline-computation time was 108s, 318s and 569s for the UR5, Sawyer and UR5 with mobile base, respectively. It should be noted, however, that this could be substantially sped up through various optimisations such as parallelisation of the candidate $\epsilon$-GHA map search (Alg. 1 lines 8-14) and IK solution computation.

In this section we discuss the experimental results and provide an analysis on HAP’s performance compared to the baselines. A key benefit of HAP is its ability to effectively reason about low-level motion during task sequencing. A notable supporting observation of this is that HAP’s task sequences tended to minimise movement between bookshelf rows. This is in contrast to the baselines which ignore the cost of avoiding collision with the shelves and objects, leading to frequent switching between rows and resulting in inefficient trajectory sequences and higher planning failure rates. See Appendix -C for supplementary videos with example experiments.

More specifically to RoboTSP, its decoupled sequencing approach fails to reason about the non-linear relationship between task space and configuration space; that is, short paths in task space are not necessarily so in configuration space. This is evidenced in the execution time plot for the UR5 experiments (Fig. 9b) where it is expected that average execution times should decrease as the number of tasks increase due to the increasing spatial density of the tasks; however, the opposite is observed for RoboTSP. An example consequential behaviour is depicted in Fig. 11 where unnecessarily long trajectories are executed.

Cluster-RTSP performs overall better than RoboTSP in the UR5 experiments. Compared to HAP, however, trajectory jerks and planning time and success rates were overall worse. A contributor of this, as earlier claimed, is that outlier tasks with highly dissimilar IK solutions can have an overall negative impact on the trajectory sequence, as depicted in Fig. 12. This can be explained by the initial unique configuration assignment step which attempts to find configurations that are all close to each other.

In contrast, the $\epsilon$-GHAs computed by HAP account for nonlinearities by considering prior knowledge such as the manipulator’s kinematic structure and environment obstacles. This is exemplified in Fig. 13 where trajectories between intra-subspace tasks appear to be consistently short and smooth. Furthermore, tasks requiring a large change in configuration are grouped into the same subspace and executed together without negatively affecting other tasks in the sequence.

Enabling the base of the UR5 robot to be mobile further enhanced performance due to the greater flexibility afforded by the diverse subspace assignments (see Figs. 15-16). Emergent behaviour such as translating the base away from poses in the lower shelf in order to avoid self-collision were observed. See Appendix -C for a supplementary video demonstrating this.

While the above effects are less pronounced in the Sawyer experiments due to its kinematic redundancy (see Fig. 14), Cluster-RTSP’s performance was notably worse than RoboTSP and HAP. This can be explained by the algorithm’s inability to effectively reason over the larger number of redundant IK solutions afforded by the arm’s extra DOF. Similar findings were reported in the original paper [4].

Finally, one may suspect that HAP produces similar goal configurations to the baselines and the seed trajectories retrieved from the $\epsilon$-GHA mappings led to higher success rates. However, we elucidate this claim by showing that HAP-no-prior still performs markedly better even without informed trajectory priors.

Here we discuss limitations of our framework and the current implementation. By reducing the search space of the planner, we intentionally trade off flexibility in the range of solutions for computation speed. However, this may lead to pruning of potentially optimal solutions. One has some control over this by, for example, utilising an increasing number of redundant $\epsilon$-GHAs to cover more of the configuration space. In the limit the full solution space can be recovered.

Furthermore, the $\epsilon$-GHA objective in (4) is a combination of task-space coverage and degree-of-path preservation which is not directly related to the task sequencing objective; the total sequence cost. Thus, it is difficult to draw any guarantees on the overall path sequence optimality. One potential work around is to store a larger number of subspaces with varying coverage of the configuration space and select a subset of these based on the online scenario.

The online planning process additionally introduces suboptimality to the overall sequence cost. For example, tasks are assigned greedily to a subspace, subspace groups are executed in arbitrary order, intra-subspace sequences are subject to the suboptimality bounds of the TSP solver and the utilisation of a fixed home configuration for transitioning between subspaces. We chose these as they were straight-forward implementations and in practice performed well in our problem setting. In general, we envision many potential approaches to utilising the $\epsilon$-GHA mappings and our online planning method is just one such approach rather than a fundamental component of the framework.

Our results and limitations discussed above motivate several important avenues of future work. For example, formulating the $\epsilon$-GHA objective in a way that directly considers the task sequence objective would be an interesting extension. Regarding the current implementation, HAP consists of several hyperparameters that may be cumbersome to tune. Automatic tuning of these could help improve usability. For example, one could train a learning model to predict $\epsilon$ based on the robot model, task space and environment.

It would be interesting to explore methods that would adapt the subspaces online in response to changes in the environment, potentially using online domain adaptation techniques [44] and conditional density estimation techniques [45]. While we focus on task sequencing in this work, $\epsilon$-GHAs are potentially useful for other applications. For example, local reactive controllers such as RMPflow [46] would benefit from knowledge of which regions of the task space are approximate isomorphisms to configuration space.

Another interesting avenue of future work is to compute the maps using a continuous representation, for example a Gaussian process, while maintaining the bounded suboptimality guarantees. Additionally, applying HAP to other embodiments such as high-DOF humanoids presents interesting challenges.