Hybrid Planning for Dynamic Multimodal Stochastic Shortest Paths

TL;DR

This paper introduces DMSSPs, a new class of MDPs for complex decision-making under uncertainty, and proposes HSP, a hybrid algorithm that improves planning efficiency and solution quality in dynamic multimodal environments.

Contribution

The paper presents a novel class of MDPs called DMSSPs and a hybrid planning algorithm HSP that combines heuristic search, approximate dynamic programming, and hierarchical planning.

Findings

HSP outperforms state-of-the-art algorithms in autonomous multimodal routing.

HSP achieves higher quality solutions with efficient computation.

The approach effectively handles uncertainty and dynamic updates in complex environments.

Abstract

Sequential decision problems in applications such as manipulation in warehouses, multi-step meal preparation, and routing in autonomous vehicle networks often involve reasoning about uncertainty, planning over discrete modes as well as continuous states, and reacting to dynamic updates. To formalize such problems generally, we introduce a class of Markov Decision Processes (MDPs) called Dynamic Multimodal Stochastic Shortest Paths (DMSSPs). Much of the work in these domains solves deterministic variants, which can yield poor results when the uncertainty has downstream effects. We develop a Hybrid Stochastic Planning (HSP) algorithm, which uses domain-agnostic abstractions to efficiently unify heuristic search for planning over discrete modes, approximate dynamic programming for stochastic planning over continuous states, and hierarchical interleaved planning and execution. In the domain of autonomous multimodal routing, HSP obtains significantly higher quality solutions than a state-of-the-art Upper Confidence Trees algorithm and a two-level Receding Horizon Control algorithm.

Figures

Tables

Table 1. Table 1 : The relevant parameter values for the three variants of the UCT baseline.

Variant	depth	exploration	n_iterations	init_N
UCT1	$100$	$50.0$	$500$	$1$
UCT2	$200$	$1.0$	$500$	$1$
UCT3	$100$	$100.0$	$500$	$50$

Table 2. Table 2 : The three sets of problems (with 1000 1000 1000 episodes each) we used for evaluating the performance of our HSP approach against the two baselines. The two parameters varied are the number of active car routes and the probability of a perturbation (speedup/delay) to the ETA of a future route waypoint at each timestep. As an additional note, each route has an average of 10 10 10 waypoints; multiplying the number of routes by 10 10 10 will give an idea of the number of route vertices being considered.

Episode Set	Num. Routes	Perturbation Prob.
Set 1 (in main)	$100\text{-}1000$	$0.75$
Set 2	$500\text{-}2500$	$0.75$
Set 3	$500\text{-}2500$	$0.35$

Table 3. Table 3 : The average number of mode switches for the various algorithms on the three problem sets, where Set 1 is reproduced from the main body for completeness . There is a sharp relative increase from Set 1 to Set 2 as the number of transit routes increases significantly. For most algorithms, there is also a slight increase from Set 2 to Set 3 as the waypoint perturbation probability is lower, i.e. the context changes less dynamically, allowing slightly more connections to be made.

Algorithm	Set 1	Set 2	Set 3
$\beta 0.75$ (HSP)	$3.422$	$4.305$	$4.399$
UCT1	$0.024$	$0.082$	$0.096$
UCT2	$0.018$	$0.037$	$0.032$
UCT3	$0.015$	$0.261$	$0.289$
RHC	$3.508$	$4.664$	$4.782$

Table 4. Table 4 : The problem set used for evaluating the computation is the first set of episodes, where the average number of car routes is 100 100 100 to 1000 1000 1000 for an episode. The decision frequency of HSP (that of its local layer except after asynchronous interrupts) is comparable to that for RHC and slightly faster than UCT. The global planning layer of HSP is slightly slower than that RHC’s global layer (the cost-to-go lookup is a bit more expensive than the nominal edge weight).

Algorithm	Avg. Time ($\mathrm{ms}$)
HSP Global Layer	$100\text{-}200$
HSP Local Layer	$10\text{-}40$
UCT1	$80\text{-}90$
UCT2	$85\text{-}95$
UCT3	$90\text{-}125$
RHC Global Layer	$100\text{-}150$
RHC Local Layer	$8\text{-}30$

Equations

$$\underset{d_{1:H-1},a_{0:H-1}}{\mathrm{argmin}}$$

$$\displaystyle c_{t+1}\sim T(d_{t},c_{t},a_{t})\$$

$$\displaystyle s_{t+1}=T\left(\xi_{t},s_{t},a_{t}\right)\$$

$$s_H\in {\mathcal{S}}_G$$

$${\zeta }_t=\left(d^0,c^0\right)\stackrel{{\pi }_{d^0}}{\to }\left(d^0,c_p^0\right)\stackrel{a\in {\mathcal{A}}_{\mathcal{D}}}{\to }\left(d^1,c_e^1\right)\stackrel{{\pi }_{d^1}}{\to }\left(d^1,c_p^1\right)\stackrel{a\in {\mathcal{A}}_{\mathcal{D}}}{\to }\left(d^2,c_e^2\right)⇝\left(d_G,c_G\right)$$

$$\text{textsc}EdgeWeight(P_{1:K},⟨d_{\text{top}},c_{\text{top}}⟩\to ⟨d^{\prime },c_p⟩)={\sum }_{k=1}^{K}P(k)J_{d_{\text{top}}}^k(⟨c_{\text{top}},c_p⟩),k=1\dots K.$$

$$T_d(⟨c,c_g⟩,a,⟨c^{\prime },c_g⟩)=T(d,c,a,c^{\prime });R_d(⟨c,c_g⟩,a,⟨c^{\prime },c_g⟩)=R(d,c,a,c^{\prime })$$

$$J_d^k(\hat{s})=\underset{a\in {\mathcal{A}}_{\mathcal{C}}}{\min }{J}_d^k(\hat{s},a)=\underset{a\in {\mathcal{A}}_{\mathcal{C}}}{\min }{\sum }_{{\hat{s}}^{\prime }\in \mathcal{C}\times \mathcal{C}}{T}_d(\hat{s},a,{\hat{s}}^{\prime })\left[-R_d(\hat{s},a,{\hat{s}}^{\prime })+J_d^{k-1}({\hat{s}}^{\prime })\right]$$

$$a_t={\pi }_{d_t}(\hat{s})=\underset{a\in {\mathcal{A}}_{\mathcal{C}}}{\mathrm{argmin}}{J}_{d_t}^{0:\bar{K}}(\hat{s})=\underset{a\in {\mathcal{A}}_{\mathcal{C}}}{\mathrm{argmin}}{\sum }_{k=1}^{K}P(k)J_{d_t}^k(\hat{s},a).$$

$${\varphi }_d=\underset{\mathbf{a}\in {\mathcal{A}}_{\mathcal{C}}^K,c\in \mathcal{C}}{\max }-{\sum }_{k=1}^{K}R(d,c,{\mathbf{a}}_k),$$

$$\text{if}{\sum }_{k=1}^{K}P(k)J_d^k(\hat{s})>{\beta }_d\cdot {\sum }_{k=1}^{K}P(k){\bar{J}}_d(k)\text{then}{a}_t\leftarrow \text{ Closed-loop Pre-emption}$$

$$J_0(c_0)=0;J_0(c\ne c_0)={\varphi }_d.$$

$${\rho }_{\pi }(c)=P\{c_{t+k}=c_0\mid c_t=c,\pi \}$$

$${\hat{J}}_d^k(c\mid \pi )=\sum _{c^{\prime }\in \mathcal{C}\times \mathcal{C}}{T}_d(c,\pi (a),c^{\prime })\left[-R_d(c,\pi (a),c^{\prime })+J_d^{k-1}(c^{\prime }\mid \pi )\right],J_d^0=0\forall c\in \mathcal{C}$$

$$J_d^k(c\mid \pi )$$

$$={\hat{J}}_d^k(c\mid \pi )+[1-{\rho }_{\pi }(c)]\cdot {\varphi }_d,$$

$${\pi }^{\ast }={\mathrm{argmin}}_{\pi \in \Pi }{\hat{J}}_d^k(c\mid \pi )+[1-{\rho }_{\pi }(c)]\cdot {\varphi }_d.$$

$$\rho (c,a)=\underset{c^{\prime }\in \mathcal{C}}{\max }\frac{\Vert c^{\prime }\Vert }{\Vert c\Vert }\text{where}{T}_d(c,a,c^{\prime })>0$$

$$\forall c\in \mathcal{C},\exists a\in {\bar{\mathcal{A}}}_{\mathcal{C}}\text{such that}\rho (c,a)<1$$

$$a_{\rho }^{\ast }(c)=\underset{a\in {\bar{\mathcal{A}}}^{\rho }(c)}{\mathrm{argmax}}\rho (c,a).$$

$$\delta =\underset{c\in \mathcal{C}}{\min }\rho \left(c,a_{\rho }^{\ast }(c)\right).$$

$$\delta^{K_{d}}<\epsilon\implies\framebox{$K_{d}=\lfloor\log_{\delta}\epsilon\rfloor$}$$

$$d_t\in \{\text{textsc}Move,\text{textsc}Ride\};c_t=\left(x_t,y_t,{\dot{x}}_t,{\dot{y}}_t\right);s_t=\left(d_t,c_t\right)$$

$$a_{\mathcal{D}}\in \{\text{textsc}Board,\text{textsc}Alight\};a_{\mathcal{C}}=\left(\ddot{x},\ddot{y}\right)$$

$${\xi }_t=\left({\mathbf{p}}_t^1,{\mathbf{w}}_t^1,{\mathbf{p}}_t^2,{\mathbf{w}}_t^2,\dots {\mathbf{p}}_t^{n_t},{\mathbf{w}}_t^{n_t}\right)$$

$${\mathbf{p}}_t^i=\text{Current position of vehicle i}$$

$${\mathbf{w}}_t^i=(({\mathbf{p}}_{1,t}^i,{\tau }_{1,t}^i),({\mathbf{p}}_{2,t}^i,{\tau }_{2,t}^i),\dots )$$

$$c_{t+1}=f(c_t,a_{\mathcal{C},t}+ϵ),ϵ\sim \mathcal{N}\left(0,\mathrm{diag}({\sigma }_{\ddot{x}},{\sigma }_{\ddot{y}})\right)$$

$$T(\text{textsc}Move,c_t,a_{\mathcal{C},t})\sim f(c_t,a_{\mathcal{C},t}+ϵ)$$

$$T(\text{textsc}Ride,c_t,a_{\mathcal{C},t})\sim {\xi }_{t+1}$$

$$T({\xi }_t,(\text{textsc}Move,c_t),\text{textsc}Board)=(\text{textsc}Ride,c_{t+1})$$

$$T({\xi }_t,(\text{textsc}Ride,c_t),\text{textsc}Alight)=(\text{textsc}Move,c_{t+1})$$

$$R(s_t,a,s_{t+1})=-(\underset{\mathrm{energy}}{\underbrace{{\lambda }_d\Vert c_{t+1}-c_t{\Vert }_2+{\lambda }_h\text{mathds}{1}_h}}+\underset{\mathrm{time}}{\underbrace{1}})$$

Full text

1 Introduction

Consider the problem of a robot arm picking and arranging objects from a conveyor belt. This simple example captures several challenges of sequential decision-making for robotics: (i) the system state is a hybrid of discrete logical modes (is the robot holding an object or not) and continuous robot state values (joint angles); (ii) external information that constrains the inter-modal transitions (the object positions that define where they can be picked up) is dynamically changing; (iii) the objective is to reach a goal (a given arrangement) with minimum cumulative trajectory cost, i.e. a stochastic shortest paths problem [1], which is challenging because the cost of a solution depends on both the high-level sequence of objects grasped and the underlying motor control actions. We define a class of Markov Decision Processes with the above properties, which we call the Dynamic Multimodal Stochastic Shortest Path (DMSSP) problem. DMSSPs can represent many general decision-making problems of importance in robotics, such as task and motion planning for mobile manipulation [2, 3] and autonomous multimodal routing [4].

Work on such robotics planning problems (hybrid state space, uncertainty, online information, solution quality as the objective) have largely focused on efficiently solving deterministic variants, delegating management of uncertainty to a low-level controller and replanning when it fails. For a time-constrained setting, however, uncertainty in the dynamics may have significant downstream effects, e.g. not reaching an object in time may invalidate the plan and greatly increase cost. A framework that accounts for uncertainty in a high level plan will be more robust to anticipating and avoiding downstream errors and temporal constraint violations. We are motivated by designing such a framework while mitigating the inevitable increase in complexity due to considering uncertainty.

Existing planning algorithms for solving large MDPs, which do account for uncertainty, would encounter several hurdles with DMSSPs. For offline MDP methods based on value or policy iteration, even with hierarchical decomposition [5] and state-of-the-art function approximators [6], it is typically infeasible to generate a good policy over the entire space of external information. Online MDP methods based on stochastic tree search suffer from the large depth and branching factor for long-horizon search over continuous states and discrete modes.

Our work builds upon the idea that a principled composition of classical search-based planning, planning for MDPs, and hierarchical planning can reason over long horizons, explicitly account for the underlying uncertainty, and replan efficiently online. We develop an algorithmic framework with three broad components: a global open-loop layer that plans for a sequence of discrete modes, a local closed-loop layer that controls the agent under uncertainty through the modes, and hierarchical interleaved planning and execution to adapt to dynamic external information. We expect the resulting approach to achieve better quality solutions on DMSSPs than existing MDP methods.

Our contributions are as follows: (i) We introduce and formulate the problem of Dynamic Multimodal SSPs. (ii) We design a Hybrid Stochastic Planning algorithm for decision-making in DMSSPs, which uses a careful choice of domain-independent representations and abstractions to efficiently incorporate heuristic search, approximate dynamic programming, and hierarchical interleaved planning and execution. (iii) We demonstrate how our approach outperforms two complementary baselines—a state-of-the-art Upper Confidence Trees algorithm and a two-level Receding Horizon Control algorithm—on real-time multimodal routing problems.

Related Work Overview

We provide a brief summary of the concepts we build upon. Our formulation is based on Markov Decision Processes (MDPs) [7] and Stochastic Shortest Paths [1], an undiscounted goal-directed negative-reward MDP. Our DMSSPs model share elements with previously studied MDP models: arbitrarily modulated transition functions [8], stochastic shortest paths with online information [9], and factored hybrid-space MDPs [10]. Our HSP algorithm uses ideas from heuristic search [11, 12] and search-based planning for multi-step tasks [13, 14], approximate dynamic programming [15, 6], hierarchical planning for solving large MDPs [16, 17, 18], and interleaved planning and execution [19, 20]. A body of relevant previous work incorporates heuristic search and classical AI techniques in algorithms for solving MDPs [21, 22, 23]. Several works from the robotics planning community solve related problems using local information to inform global planning and trajectory optimization [24, 25] and explore various aspects of combined (discrete) task and (continuous) motion planning [3, 2] such as planner-agnostic abstractions [26], stochastic shortest path formulations [27, 28], and hierarchical planning and execution [29, 30]. Our optimization-based formulation is similar to that for Logic-Geometric Programming [31].

2 Dynamic Multimodal Stochastic Shortest Paths (DMSSPs)

A DMSSP is a discrete-time MDP, $\mathcal{M}=(\mathcal{S},\mathcal{A},\Xi,T,R)$. The state space is multimodal, i.e. factored as $\mathcal{S}\equiv\mathcal{D}\times\mathcal{C}$. Each $d\in\mathcal{D}=\{1,\ldots,{\absolutevalue{\mathcal{D}}}\}$ represents a discrete logical mode of the system and $\mathcal{C}$ is the continuous state space of the decision-making agent. The current system state is accordingly denoted as $s_{t}=\left(d_{t},c_{t}\right)$. The action space is $\mathcal{A}\equiv\mathcal{A}_{\mathcal{D}}\cup\mathcal{A}_{\mathcal{C}}$ where $\mathcal{A}_{\mathcal{D}}$ is the set of discrete mode-switching actions (e.g. Pickup and Putdown) and $\mathcal{A}_{\mathcal{C}}$ is the set of control actions for the agent. The context space $\Xi$ is dynamic. At each time-step, the agent observes a discrete set of contexts $\{\xi_{t},\xi_{t+1\mid t},\xi_{t+2\mid t},\ldots,\xi_{t+K\mid t}\}$. Here, $\xi_{t}$ is the current context and $\xi_{t+k\mid t}$ is the estimated context $k$ time-steps in the future. In general, $\xi_{t+k\mid t}\neq\xi_{t+k\mid t+1}$. The context is assumed to be generated by an exogenous process whose evolution equation, i.e. $\xi_{t+1}=g(\xi_{t})$, is complex and depends on a number of unobserved variables. The current context set is compactly denoted as $\xi_{t:t+K}$. The set cardinality $K$ is some domain-dependent prediction horizon, and the dimensionality of context is time-varying, i.e. $|\xi_{t}|$ may differ from $|\xi_{t+1\mid t}|$. In our example, the context is the current and estimated future positions of all moving objects on the belt.

The transition function $T$ can be factored as follows:

[TABLE]

The factored form efficiently encodes stochastic intra-modal physical dynamics in eq. 1a and context-dependent deterministic mode-switching rules in eq. 1b. We treat the discrete mode-switching as single time-step actions. The grounded values of logical predicates (for mode-switching) are represented as states [27], i.e. instances of $\mathcal{S}$ and not just $\mathcal{D}$ in eq. 1b. The current context $\xi_{t}\in\Xi$ restricts the actual grounding of the logical predicates in the preconditions and effects that define the mode switch rules. Precondition states need to be reached for a mode switch to be feasible and effect states are obtained after the mode switch is completed. In our running example, whether the end-effector can satisfy the preconditions to grasp a moving object depends on the end-effector (state space) and the object (context). As in typical stochastic shortest paths [1], the problem is episodic and undiscounted. Also, the reward function is non-positive, so we will also use ‘cost’ to represent negative rewards. It is factored in terms of the intra-modal costs, i.e. $R:\mathcal{D}\times\mathcal{C}\times\mathcal{A}_{\mathcal{C}}\times\mathcal{C}\to\mathbb{R}^{-}$. The decision diagram of a DMSSP is depicted in fig. 1(a).

We explicitly write out a DMSSP as an online stochastic optimal control problem. Given a start state $s_{0}=(d_{0},c_{0})$ and goal space $\mathcal{S}_{G}=(d_{G},\mathcal{C}_{G}\subseteq\mathcal{C})$, the overall problem is:

[TABLE]

where $H$ is the number of steps taken. From fig. 1(b), a DMSSP planning algorithm has to choose a sequence of modes from start mode $d_{0}$ to goal mode $d_{G}$, and also control the agent to traverse through the corresponding subspaces. Even the one-shot deterministic problem is a discrete-time Logic-Geometric Program [31], for which a globally optimal solution is infeasible. Theoretical optimality analysis requires several assumptions and is not our focus (additional comments in appendix A). We are interested in high quality solutions that are efficient.

3 Hybrid Stochastic Planning (HSP)

The challenges with DMSSPs preclude directly applying existing techniques for large MDPs, as mentioned in section 1. We develop a Hybrid Stochastic Planning framework that is particularly suited for DMSSPs. A global layer computes a sequence of mode switches and corresponding transition states; planning open-loop enables efficiently searching over long-horizon mode sequences. A local layer executes actions that control the agent within each mode; closed-loop planning provides some robustness to uncertainty. Additional hierarchical logic for interleaving planning and execution reacts to the dynamically changing context at both global and local levels.  Figure 2 depicts a schematic of the HSP structure. The concepts we rely on have been studied extensively. However, what we particularly contribute are the design choices in the overall framework (discussed subsequently); heuristic search over mode sequences, decomposition to smaller MDPs, using cost-to-go functions, i.e. negative value functions as edge weight surrogates, pre-empting the local controller. These choices unite techniques from classical search-based planning, planning for MDPs, and hierarchical planning in a principled manner. Since DMSSP is a class of problems, we present a general algorithmic framework that can accommodate a variety of modular techniques and subroutines.

3.1 Global Open-loop Layer

This layer repeatedly computes a high-level plan from the current state $s_{t}=(d_{t},c_{t})$ to a goal state $(d_{G},c_{G}\in\mathcal{C}_{G})$. The current high-level plan (computed at current time $t$) is denoted as follows:

[TABLE]

where $(d^{0},c^{0})=(d_{t},c_{t})$ and the sequence of chosen modes is $d^{0}\rightarrow d^{1}\rightarrow\ldots\rightarrow d_{G}$. For each planned mode switch $d^{n}\rightarrow d^{n+1}$, the sampled precondition is $c_{p}^{n}$ and effect is $c_{e}^{n+1}$. For each mode, $\pi_{d^{n}}$ is a closed-loop policy that controls the agent through the region from the effect after the previous switch $c_{e}^{n}$ to the precondition before the next switch $c_{p}^{n}$. GlobalPlan in Algorithm 1 outlines the global layer, which runs a heuristic search [11, 13] from the current state $s_{t}$ to compute a sequence of modes and transition states towards the goal space $\mathcal{S}_{G}$. We briefly discuss the subroutines.

For NextValidModes, we use the mode-switch rules and the current context set to identify possible next modes from $d_{\text{top}}$. The specific implementation is domain-dependent. In our example, the robot can only execute Pickup on objects that will be reachable based on the current context set. The SampleTransition method generates grounded precondition and effect states for the proposed mode switch given the context set, using some domain-specific sampling procedure (the number of samples $N$ is a parameter). It also generates $P_{1:K}$, a probability distribution over future time-steps at which the context $\xi_{t+k}$ will satisfy the mode transition preconditions. In our example, we could use work on sampling end-effector poses for grasping a given object [32] at various steps along the object’s expected future trajectory (defined by the context).

For EdgeWeight, consider its role in the overall solution. Each call to GlobalPlan returns a sequence of modes and transition states. Subsequent execution by the local layer can only optimize the trajectory locally. The edge weight function should therefore be a good surrogate of the expected cost to have any hope of a good overall solution. The cumulative intra-modal dynamics cost depends on the actual trajectory that will eventually be executed by the agent in that mode. Therefore, as a surrogate, we use the cost-to-go function of the policy for a local MDP corresponding to the mode, where the local MDP state is encoded in the edge (details explained in section 3.2):

[TABLE]

Here, $J^{0:K}_{d_{\text{top}}}$ is the horizon-dependent cost-to-go function for the local MDP of mode $d_{\text{top}}$. We have presented the most general search formulation for the global layer. It can use many of the speedup techniques for heuristic search-based planning [14, 33] to improve efficiency. The heuristic function ($\mathrm{heur}$) in general is defined on $\mathcal{D}\times\mathcal{D}$; additional comments on the heuristic are in appendix B.

3.2 Local Closed-loop Layer

The local closed-loop layer in HSP controls the agent in its current mode up to the chosen transition state for the next switch. The layer is ‘local’ because it is only provided information about the currently executing step of the current global plan $\zeta_{t}$, i.e. $(d^{0},c^{0})\rightarrow(d^{0},c_{p}^{0})$, where $d^{0}$ is the current mode of the agent. For each mode $d$, we define a local target-directed MDP $(\mathcal{C}\times\mathcal{C},\mathcal{A}_{\mathcal{C}},T_{d},R_{d})$, where the first argument of the state is the current position, and the second is the target. The control space is $\mathcal{A}_{\mathcal{C}}$. The $T_{d}$ and $R_{d}$ functions are derived from $T$ and $R$, i.e.

[TABLE]

The Cartesian self-product $\mathcal{C}\times\mathcal{C}$ is required in general because the policy must be able to control the agent between any two states in $\mathcal{C}$. However, for many spaces, a state $\langle c_{1},c_{2}\rangle$ can be encoded with relative state $c_{1}\circ c_{2}\in\mathcal{C}$, where $\circ$ is a difference operator, and the target is always the ‘origin’ $c\circ c$, i.e. the zero state. Many modes typically share dynamics, so the same policy can be reused [17]. For example, in our setting, the robot dynamics effectively depend only on if an object is currently grasped or not, which can be encoded with two modal dynamics functions.

**Finite-Horizon Value Iteration

**The local closed-loop policy has a dual role, controlling the agent with low cost within the mode to the transition state $c_{p}$ chosen by the global layer and satisfying the temporal constraints of the context for the next mode switch ($P_{1:K}$). If it is too slow, the mode transition may fail, affecting the overall solution. On the other hand, controlling the agent as quickly as possible may be highly sub-optimal. To model this tradeoff, we use finite-horizon value iteration to obtain a horizon-dependent policy [15]. The finite-horizon value iteration requires a horizon limit (we use the context horizon $K$) and a terminal cost $J^{0}_{d}$. For all local MDP states $\hat{s}\in\mathcal{C}\times\mathcal{C}$, we compute for $k=1\ldots K$,

[TABLE]

where $J_{d}^{k}$ represents both cost-to-go and action-cost-to-go (Q function), with overloading (negative reward i.e. positive cost). The full cost-to-go function, compactly denoted as $J^{0:K}_{d}$, can be used in eq. 4. As in the global layer, we require only a general framework for value iteration; any local or global approximation scheme [6] and other approximate dynamic programming [15] techniques could be used. The regional closed-loop policy $\pi_{d}$, invoked during execution online, is based on $J^{0:K}_{d}$ obtained offline. For the current DMSSP state $(d_{t},c_{t})$, context horizon distribution $P_{1:K}$, and transition state $c_{p}$ (provided by the global layer), the control action chosen locally for $\hat{s}=\langle c_{t},c_{p}\rangle$ is

[TABLE]

**Terminal Pseudo-Cost

**To incentivize the closed-loop policy to reach the target, we need a terminal penalty for states where the target is not reached at horizon 0, i.e. $\{\langle c,c_{g}\rangle\in\mathcal{C}\times\mathcal{C}\mid\lVert c-c_{g}\rVert>\epsilon_{d}\}$, for some domain-dependent distance metric. We need to be careful while choosing the terminal penalty or pseudo-cost due to the sub-optimality of hierarchical MDP planning [5, 18]. The penalty is not in the true cost function in eq. 2, so the higher it is set, the poorer is $J^{0:K}_{d}$ as a surrogate of the true cost, and the more (locally) sub-optimal is $\pi_{d}$ as a controller. An insufficiently high penalty, on the other hand, may lead to $\pi_{d}$ choosing lower-cost actions at the risk of being unable to reach the target within the context horizon. Consequently, the attempted mode switch may fail, forcing HSP to recompute a different mode sequence, leading to much poorer solutions (downstream effects of uncertainty).

We set the penalty $\phi_{d}$ as the maximum cost of any $K$-length action sequence within the mode, i.e.

[TABLE]

The pseudo-cost is the smallest penalty that prioritizes the mode sequence. The $\phi_{d}$ value can be computed offline and then used in the finite-horizon value iteration. Further comments on the pseudo-cost and horizon limit are made in appendices A and C). LocalPreprocessing in algorithm 1 outlines the local layer; the policies obtained from pre-processing are used online.

An implicit assumption of ours is that the MDP for the agent dynamics can be solved reasonably. This assumption is not always valid, e.g. a complex underactuated system or an articulated manipulator. However, for many practical systems, framing and solving the control problem with MDPs has been successful [34], and finite-horizon versions of those controllers could also be used here.

3.3 Hierarchical Interleaved Planning and Execution

Interleaving planning and execution is an important property for real-time decision-making. Our HSP framework uses the simplification approach to interleaving [19]. The global layer simplifies the underlying intra-modal control problem by determinization and planning over multiple timesteps, and computes a solution in this simplified space. The local layer executes the plan provided by the global layer. We discuss two aspects of our interleaving, each occurring at one of the levels, either global replanning or local pre-emption. The full HSP framework is outlined in algorithm 2.

**Global Replanning

**HSP uses a combination of event-driven and periodic replanning [35]. The two events that trigger replans asynchronously are closed-loop pre-emption and failed mode-switch attempts. The current global plan is then invalidated and a new plan must be generated before execution can resume. With periodic replanning, the global layer computes a new global plan from the current state synchronously in the background, while the local layer executes the current plan. The duration of the replanning period is domain-dependent; but a domain-agnostic strategy is to replan immediately after the previous plan has finished computing. This updates the local layer’s next target every $\Delta T$ time-steps, where $\Delta T>1$ is the ratio between closed-loop and open-loop frequencies. Each time a global plan has been recomputed by the global layer, the next target of the local policy $\pi_{d}$ is immediately updated to the (potentially) new transition state for the next switch (fig. 3 illustrates this).

**Local Pre-emption

**Updates to $\xi_{t:t+K}$ can make it difficult to reach the chosen mode transition state in time (e.g. a target object suddenly speeds up). The periodic global replanning does account for this. However, the latency is $\Delta T$ time-steps. We have additional closed-loop pre-emption logic to reason about the next chosen mode switch at the higher frequency of the local layer. For each local MDP, we compute and store (during pre-processing) the worst cost-to-go for any state from each horizon value, i.e. $\bar{J}_{d}(k)=\max_{\hat{s},a}J^{k}_{d}(\hat{s},a),\ k=1\ldots K$. During the online execution of the local layer, if the agent is at a state from where reaching the goal in the remaining horizon is sufficiently unlikely, i.e.

[TABLE]

where $\beta_{d}\in[0,1]$ is a risk parameter. The lower $\beta_{d}$ is set, the lower the risk we are willing to take that the agent can reach the next transition state in time. In addition to a higher frequency for reasoning about pre-emption, this logic provides a modifiable risk aversion.

4 Experiments: Multimodal Routing Domain

We use a different domain than our running example; the recently introduced Dynamic Real-time Multimodal Routing (DREAMR) problem [4]. We omit an elaborate description of the domain (see appendix E for more). The DREAMR problem requires planning and executing routes under uncertainty for an autonomous agent that can use multiple modes of transportation in a dynamic transit vehicle network. There are two discrete modes in the problem, Move for when the agent moves by itself and Ride for when the agent uses transport. The continuous state is the agent’s position and velocity. The mode-switching actions are Board and Alight, which switch Move to Ride and vice versa respectively. The noisy agent control actions are for acceleration in each direction. The agent is penalized for energy expended due to movement and waiting in place, and total elapsed time. The transit vehicle routes are the contextual information; at any time, the current context set comprises the current position and estimated future route (as a sequence of waypoints) of each active vehicle. The estimated time of arrival (ETA) for each subsequent waypoint is subjected with some probability to a bounded two-sided deviation at each timestep. Mode switches can only be made at transit waypoints. Additionally, the agent can only Board a vehicle if it is sufficiently close to the waypoint at the same time as the vehicle, and sufficiently slow. The dimensionality of the context space increases with the number of active route waypoints, and so is highly dynamic and very large (in the thousands). Though there are technically only two modes, in practice, the number of valid mode sequences to the goal is exponential in the number of transit vehicle routes.

4.1 Baselines: Upper Confidence Trees and Receding Horizon Control

We use two complementary baselines to evaluate the benefit of our solution. First, a (domain-specific) two-level Receding Horizon Control (RHC) method which repeatedly solves a deterministic problem. It uses graph search for planning routes, and non-linear receding horizon control trajectory optimization for executing them. The second is based on Upper Confidence Trees (UCT) [36], a general online MDP planning algorithm, with two enhancements: (i) techniques from PROST [37], a state-of-the-art UCT-based probabilistic planner [38]; (ii) double progressive widening [39] which artificially limits the branching factor and is more suitable for a continuous state space. We additionally assist the UCT baseline: (i) the value function estimates and Q-value initializations are informed by $J^{0:K}_{d}$ (ii) the tree depth is set to the same horizon limit $K$ as the HSP local layer (iii) many trials are run to compute good estimates for each action.

4.2 Results

We used the POMDPs.jl framework [40] in Julia (additional details in appendix D and attached code111The Julia code is at https://github.com/sisl/CMSSPs). For HSP, we chose three values for the $\beta\in[0,1]$ risk parameter for attempting time-constrained risk parameters: $0.55$, $0.75$ and $0.95$ (for lower values, it is overly risk-averse and rejects most transit connections). For UCT, we chose three different combinations of the number of virtual rollouts at a new tree node and the UCB exploration constant. The RHC does not have an important tuning parameter. We used the same large-scale problem scenarios as the original DREAMR paper (see appendix E for more details and parameter values). We evaluated each algorithm instantiation on $1000$ simulated episodes.  Figure 5 shows the average cumulative trajectory cost; both UCT and RHC produce poorer quality trajectories than HSP.  Figure 5 displays the average number of executed mode switches. Appendix E.3 has average computation times, and performance results on two additional sets of episodes with varying problem parameters. HSP has low sensitivity to values of $\beta>0.5$ (at least on this domain), which is a useful property. Multiple factors explain the relative performance gap between HSP and the baselines. For UCT, compared to HSP, it uses far fewer modes of transportation. For most episodes, UCT controls the agent directly to the goal without any mode switches. As a result, it incurs a far higher overall cost for energy expended due to movement. In general, a tree-based method such as UCT requires a thorough search over the actions with very deep lookahead to even possibly consider taking transportation, as we identified earlier. On the other hand, RHC uses about the same number of mode switches as HSP, due to the similar lookahead in its global layer enabled by deterministic planning. However, it uses a nominal edge weight for the graph search, which is a poorer surrogate than the cost-to-go function, and makes the choice of mode switches poorer. It also uses receding horizon control rather than closed-loop control at the local level, which makes it more sub-optimal locally, especially for making time-constrained connections.

5 Discussion

We introduced and formulated the Dynamic Multimodal Stochastic Shortest Path problem for representing sequential decision-making problems for complex robotics settings. Our Hybrid Stochastic Planning framework, through our choice of abstractions, is a principled way of incorporating techniques from heuristic search, approximate MDP planning, and interleaving planning and execution. HSP’s performance on the real-time autonomous routing domain against the complementary baselines highlights our general motivation. By explicitly using online long-horizon planning and accounting for the underlying uncertainty and its downstream effects, we can achieve good quality solutions. Our key limitations are the assumptions on the various subroutines and components, e.g. the mode transition states can be sampled efficiently from the context, the MDP for the agent dynamics can be solved, the value function lookup is fast, and so on. However, as we mentioned, there are still several domains of interest where these assumptions are quite reasonable and have been used effectively, and our work would be applicable in all of them. Future research includes more detailed theoretical analyses with problem assumptions, using an online stochastic planner at the local layer to overcome the need for an offline phase, and empirical results on other problem domains.

Appendix A Comments on HSP Optimality

We briefly stated in section 2 how the properties of DMSSPs, namely the online context updates and the hybrid (discrete and continuous) state space, make any useful theoretical analysis very difficult. We provide additional comments and justification for that here and analyze quantitatively the role of the terminal pseudo-cost from eq. 8. This section is intended more to highlight the issues with optimality analysis of HSP rather than prove any particular results (which would require several modeling assumptions and is not the scope of this work).

In general, analysis of an online optimization algorithm is done with respect to the best solution in hindsight, but even this comparison typically assumes a specific functional form for the online information. However, for DMSSPs the form of the context is entirely domain-dependent (e.g. route information in DREAMR or object trajectories in dynamic TAMP). Therefore, for our subsequent discussion we will assume full observability of the true context at all future timesteps.

A.1 Global Optimality

The cost of a solution depends on the discrete sequence of mode switches chosen as well as the underlying control actions. In a deterministic setting alone, for a general non-linear cost function (which is the case for DMSSPs), this is a Mixed Integer Program (MIP), which is known NP-Complete. The presence of uncertainty in the control outcome makes this an even more difficult Stochastic Mixed Integer Program (SMIP), which is out of the scope of this discussion. Practical solutions for MIPs use heuristic methods based on combinatorial techniques like tabu search, hill climbing, simulated annealing, and others.

In our case, HSP’s global open-loop layer uses an anytime search method parameterized by the number of samples $N$ for each considered mode switch (see algorithm 1). A higher value of $N$, i.e. a greater amount of computation time devoted to the open-loop planning will yield better quality solutions. Of course, the quality of these solutions is in terms of the edge weight, which is a surrogate objective for the true cost that depends on the actual executed trajectory. Therefore, even in the asymptotic case, i.e. as $N\rightarrow\infty$, it appears that no guarantee of global optimality can be made. We use the cost-to-go function of the local MDP state encoded in the edge as a good surrogate of the expected cost.

A.2 Local Optimality

In the multimodal setting of DMSSP, local optimality refers to optimality within the chosen mode sequence, i.e. whether the HSP solution has minimum expected cost out of all the solutions constrained to follow that mode sequence. Due to the continuous component of the state space for each mode, the expected cost-to-go within a region depends on the approximation error $\epsilon_{d}$ of the value iteration method used to obtain the local policy $\pi_{d}$. Additionally, the global open-loop planning would have to cover the space of all possible sampled pre-conditions for each mode switch. Therefore, we will consider the case when the agent’s component of the state space, $\mathcal{C}$ is discrete rather than continuous, and when the global layer considers every possible discrete precondition during planning. As $N\rightarrow\infty$ and $\epsilon_{d}\rightarrow 0\ \forall d\in\mathcal{D}$, this discrete case performance will be emulated.

N.B. The following discussion relies heavily on section 7.2 of Bertsekas [15], which discusses Stochastic Shortest Path Problems for the discrete state space case.

For simplicity, we assume the following properties for the agent state space $\mathcal{C}$ and every local intra-modal MDP:

• •

$\mathcal{C}=\{c_{0},\ldots,c_{n}\}$ is a set of discrete states.

• •

There is a difference operator $\circ$ such that $c_{1}\circ c_{2}\in\mathcal{C}$. Furthermore, $c\circ c=c_{0}\ \forall c\in\mathcal{C}$. We can transform the cost function as $R_{d}(\langle c_{1},c_{2}\rangle,a,\langle c^{\prime}_{1},c_{2}\rangle)=R_{d}(c_{1}\circ c_{2},a,c^{\prime}_{1}\circ c_{2})$

• •

The state $c_{0}$ is the cost-free absorbing state. We also refer to this as the ‘zero’ state or origin, for obvious reasons.

• •

There is at least one proper policy ([15], cf. Assumption 7.2.1 footnote), i.e. a stationary policy which has non-zero probability of reaching the zero state $c_{0}$ after some number of stages $m$, regardless of the initial state. This assumption is actually quite weak in practice.

Given these assumptions, the HSP local layer sub-task of reaching state $c_{e}$ from a start state $c_{p}$ is equivalent to the classical stochastic shortest paths problem of reaching the zero state $c_{0}$ from $(c_{p}\circ c_{e})$ with minimum expected cost. DMSSPs additionally have a finite-horizon setting for $d$ because of the temporal constraints induced on mode transitions by the context set. For the fully observable context set, the distribution over future timesteps $P_{1:K}$ collapses to an exact time horizon, say $k$, within which the agent must reach the terminal state, in order to successfully make the mode switch. By Proposition 7.2.1 (a) of Bertsekas [15], for the finite-horizon case, the value iteration algorithm of eq. 6 yields the optimal stage-wise, i.e. horizon-dependent cost $J^{0:K}_{d}$ from every start state, where the terminal cost is given by $J^{0}$.

Given the start state $c_{e}$ in the current mode, if the global layer samples every possible discrete precondition $c_{p}$ for the next mode-switch out of the current region, then every possible relative start state $c_{p}\circ c_{e}$ would be considered, and $c^{*}_{p}=\mathrm{argmin}_{c_{p}\in\mathcal{C}}J^{k}_{d}(c_{p}\circ c_{e})$ would be chosen, where $k$ is the known time interval for the future context permitting the mode switch.

Therefore, at least in the discrete case, our representation of the local layer’s local MDP allows us to inherit the optimality properties of SSP problems. However, due to the finite-horizon setting, we are only optimal with respect to the terminal cost $J_{0}$. The terminal cost issue illustrates the potential conflict at the local layer between reaching the target with a low cost and reaching the target in time, and is a caveat to any local optimality analysis of DMSSPs.

For a given mode switch chosen by the global layer, the sequence of actions with minimum expected cost may not have the lowest probability of reaching the chosen transition point, i.e. reaching the zero relative state in time, which would lead to globally poorer overall trajectories. In a deterministic setting, we could have constrained sets of control actions guaranteed to reach the zero state. However, for the stochastic setting of DMSSPs, we must instead consider the probability of the action sequence to reach the zero state. We attempt to balance this with the terminal pseudo-cost of eq. 8, which we analyze further here.

Terminal Pseudo-Cost Analysis

We are considering the finite-horizon value iteration of eq. 6 with the terminal cost of eq. 8. For the discrete state space case that we are analyzing, the terminal (at horizon 0) relative state corresponding to a target being reached is the ‘zero’ state $c_{0}$. All other states represent the target not being reached and are assigned the terminal penalty from eq. 8. Therefore, the terminal cost function $J^{d}_{0}$ is defined as follows:

[TABLE]

Denote (as we did earlier) the policy derived from the optimal cost-to-go function $J^{0:K}_{d}$ as $\pi^{0:K}_{d}$ (denoted $\pi$ hereafter compactness). For any given relative state $c\in\mathcal{C}$, define the probability of the ‘zero’ state being reached from it by following $\pi$ after $k$ steps as

[TABLE]

for the current time-step $t$. Furthermore, denote the k-stage expected cost-to-go for $\pi$ with terminal cost [math] for all states as

[TABLE]

where, once again, negative reward is used to imply positive cost (the reward function is non-positive). Using eqs. 10 to 12, we can express the general k-stage cost-to-go for $\pi$ as

[TABLE]

where the left-hand term is the expected cost due to the trajectory and the right-hand term is the penalty weighted by the probability of failure to reach the target in time (i.e. in $k$ steps). If we express our finite horizon value iteration as a finite horizon policy iteration (using the fact that a policy can be extracted from a value function), then our corresponding policy search is

[TABLE]

From eq. 14, the two terms of interest for the policy iteration are the expected cost of the policy $\hat{J}_{d}^{k}(c\mid\pi)$ and the probability of failure to reach the target $[1-\rho_{\pi}(c)]$. The above analysis explicitly shows how $\phi_{d}$ is a scaling factor that balances these two terms. By setting it to the quantity in eq. 8, we are effectively prioritizing policies that reach the target in time over ones that do so with lower expected cost within the region. Ultimately, this choice of a penalty term is still a heuristic.

Appendix B Global Layer Heuristic

The heuristic function $\mathrm{heur}$ in algorithm 1 is used by the global layer planning to guide the search over good mode sequences and (hopefully) make it more efficient than searching over all possible mode sequences, which may be unacceptably expensive. Heuristic functions in search [12] are usually goal-directed (usually an easy-to-compute estimate of the cost to reach the goal from the current state). They also usually operate on points, not spaces. Therefore, a heuristic function of the form $\mathrm{heur}(s_{t},\mathcal{S}_{G})$ where the second argument is a goal space will not necessarily be well-defined. Of course, domain-specific heuristics may be able to work with a goal space (for instance by sampling a goal state), but we cannot make a general comment on that.

Therefore, in algorithm 1, the heuristic operates only on the set of modes, i.e. $\mathrm{heur}:\mathcal{D}\times\mathcal{D}\rightarrow\mathbb{R}^{+}$, and for a given state $s_{t}=(d_{t},c_{t})$, the heuristic value is $\mathrm{heur}(d_{t},d_{G})$ where $d_{G}$ is the goal mode, i.e. the mode of the goal space. There has been a long line of work on domain-agnostic heuristics for classical planning [12, 14, 33], many of which can be utilized here while searching over the logical modes. As we mentioned in section A.1, global optimality guarantees cannot be made for DMSSPs, so the heuristic functions need not be admissible, i.e. under-estimates of the cost to reach the goal. In any case, due to the uncertainty in the outcome of the underlying control actions, it can in general be challenging to devise useful admissible heuristics in stochastic settings. A potentially useful non-admissible heuristic could be based on a worst-case traversal cost within modes.

Appendix C Horizon Limit Selection

In the main text, for simplicity, we assumed that the horizon limit used is the same as the context horizon limit $K$, which is an appropriate choice if that value is known for a particular domain. In this section we discuss some general issues about the parameter and one possible choice that does not require knowledge of the context horizon limit. For the subsequent discussion, we denote the local layer’s horizon limit parameter as $K_{d}$ to distinguish it from the context horizon $K$. The choice of horizon limit parameter for the local layer influences the lookahead of the local MDP policies and the amount of memory storage required for the value function (increases linearly with $K_{d}$).

A very small limit will make the local layer more sensitive to the probability distribution over context horizon $P_{1:K}$ for the next mode switch. In our conveyor belt example, suppose the next planned mode switch is to pick up a box at a point on its future expected trajectory. If $K_{d}\ll K$, then for most future points on the box’s trajectory, we will have $\sum_{k>K_{d}}^{K}P(k)\sim 1$, i.e. the bulk of the probability mass is on future horizon values greater than the local layer horizon limit. The local layer cannot then choose a useful control action from its cost-to-go function $J_{d}^{0:K_{d}}$ in eq. 7 (the $K$ from eq. 7 is $K_{d}$ in this case). It has to wait until $\sum_{k\leq K_{d}}^{K}P(k)$ is sufficiently greater than [math], but that restricts the total reaction time of the local layer. Thus, it can only choose transition points up to a few timesteps ahead, which reduces the robustness to downstream effects, i.e. missing time-constrained mode switches.

A very large horizon value will accommodate the context horizon but will make the value function $J^{0:K}_{d}$ expensive to compute and store. Subsequently, we propose a domain-agnostic strategy that does not require knowledge of the context horizon limit. We make the same assumptions on the local intra-modal MDP as in section A.2. Also denote $\lVert c\rVert$, i.e. the norm of the relative state as the distance to the origin; the target is reached when $\lVert c\rVert\leq\epsilon$ for some $\epsilon$. The intuitive idea is to set $K_{d}$ at least high enough to allow any relative state to reach the origin with some set of control actions.

Let $\bar{\mathcal{A}}_{\mathcal{C}}=\mathcal{A}_{\mathcal{C}}\setminus\text{no-op}$ be the set of all control actions excluding the no-op action. We define the worst-case progress of an action $a\in\bar{\mathcal{A}}_{\mathcal{C}}$ as

[TABLE]

where $\rho(c,a)<1$ indicates that the action makes progress towards the origin, as the distance to the origin after the action has reduced. for all possible next states. Another assumption we make is that progress towards the origin can be made from every state, i.e.

[TABLE]

and the set of all such progress actions for a given relative state is denoted $\bar{\mathcal{A}}^{\rho}(c)$. Furthermore, define the most progressive action for a relative state as the one with maximum worst-case progress, i.e.

[TABLE]

Finally, let the minimum progress from any state in one step be denoted as

[TABLE]

By the assumption in eq. 16 and by eq. 17, we know that $\delta<1$. From any relative state, there is always an action that can reduce the distance to the origin by a fraction of at least $\delta$. Therefore, a suitable choice for $K_{d}$ satisfies

[TABLE]

where $\epsilon$ is the domain-dependent threshold parameter.

The above analysis of eqs. 15 to 19 does assume that the various involved quantities are computable for a continuous space $\mathcal{C}$. For simple dynamical systems it may be possible to do so analytically without sampling states from the space, otherwise, a sampling scheme can be used during preprocessing to generate an exhaustive set of samples and then we can directly apply the equations to compute a horizon limit parameter which is sufficient for most of the space.

Appendix D Implementation Details

N.B. Due to legacy naming reasons, the attached code uses the acronyms ‘CMSSP’ and ‘HHPC’ instead of ‘DMSSP’ and ‘HSP’ respectively.

As we mentioned earlier, our implementation is in the Julia programming language. We also rely heavily on the POMDPs.jl [40] framework for modeling and solving Markov Decision Processes. Our implementation broadly consists of: (i) a domain-agnostic component which defines the DMSSP problem model and an interface for defining the various components, and the general HSP solution framework; (ii) a domain-specific component which instantiates the various DMSSP components (discrete modes and continuous state space, transition, reward, context) and the other functions required by the HSP solution framework (NextValidModes, SampleTransition, and so on). We briefly describe the domain-agnostic component here and the domain-specific component (for multimodal routing) in section E.1. We provide more elaborate technical details in the README of the attached code.

The DMSSP problem formulation is specified in the src/models/ folder. It is primarily an interface that is parameterized by the domain-specific datatypes for the state and action spaces and the context set, that the domain designer has to provide. The HSP solution framework is implemented in the src/hhpc/ folder. As for the algorithms in the paper, it is divided into a global_layer, a local_layer, and the full hhpc_framework.

The global_layer implements the GlobalPlan procedure of algorithm 1 using a modified A* search algorithm implementation [11]. The modification is for implicit graphs, i.e. where the edges are not specified before the search is called, but rather generated on-the-fly by a successor function when a node is expanded. An implicit search allows greater flexibility, especially in domains with many discrete modes; an efficient NextValidModes subroutine can be used to only generate the modes actually reachable from the current mode, rather than explicitly enumerating all of them apriori. The local_layer implements the LocalPreprocessing procedure of algorithm 1, i.e. the terminal cost computation and finite-horizon approximate value iteration. We use multilinear grid interpolation [6] of the value function over the continuous state space, however, an alternate implementation with a different approximator could also be done here.

The hhpc_framework defines the top-level behavior of algorithm 2. It assumes access to a discrete-time domain-specific DMSSP simulator. At every time-step, it observes the current state and context set, makes a decision based on the HSP framework, and outputs an action to the simulator. It also implements the hierarchical interleaving of planning at the global layer and execution at the local layer.

Appendix E Further Experimental Details

The experimental domain used is the recently introduced dynamic real-time multimodal routing (DREAMR) problem [4]. An agent (for example a drone) has to be controlled from a start to a goal location. There is a network of transit vehicle routes that the agent has access to, and it may use transit vehicles as temporary modes of transportation along segments of the routes, in addition to moving on its own to the destination. The objective is to reach the goal location while incurring as low a trajectory cost as possible, where the agent is penalized for energy expended due to distance traversed and hovering in place, and total elapsed time.

In this section, we provide a more elaborate description of the experiments we ran to evaluate HSP against the two baseline methods of two-level Receding Horizon Control (RCH) and Upper Confidence Trees (UCT) with enhancements from PROST and Double Progressive Widening. Several implementation aspects pertaining to the problem domain of Dynamic Real-time Multimodal Routing (DREAMR) were obtained from the open-source Julia repository DreamrHHP.jl (here) of the original paper. All the experiments were run on Linux with $16\text{\,}\mathrm{GiB}$ RAM and a $6$-core $3.7\text{\,}\mathrm{GHz}$ CPU.

E.1 Domain-Specific Implementation

On the DMSSP problem formulation side, the state, action, transition, and reward functions are all obtained directly from the DREAMR paper. The context is the current position and estimated remaining route (as a sequence of time-stamped waypoints) for all currently active transit vehicles. The information is summarized below:

[TABLE]

where $\mathds{1}_{h}=\mathds{1}[\lVert(\dot{x_{t}},\dot{y_{t}})\rVert<\epsilon]$ indicates hovering in place, $\lambda_{d}$ and $\lambda_{h}$ are scaling parameters for how important the distance traversed and hovering are with respect to each other and with respect to each unit of elapsed time.

On the HSP solution side, there are only a few details worth mentioning. For simplicity, the transition between modes is constrained to only happen at the transit vehicle route waypoints. For Move to Ride, the agent can Board at potentially any of the future route waypoints of all the currently active cars, while for the converse, the agent can Alight at any of the future route waypoints of the transit vehicle it is currently on. We make this simplification because, at least for this work, we are not interested in tuning the sampling parameter $N$ and evaluating its effect on the performance of HSP. In practice, we can certainly increase resolution by sampling transition points in between the pre-decided route waypoints. In any case, because the ETA at the route waypoints is perturbed with high probability at each timestep, we simulate dynamically changing contextual information.

The agent’s $(x,y)$ position is assumed to be bounded on a $1\times 1$ unit grid (that can be arbitrarily scaled to represent a real-world grid). Since the agent’s control is holonomic, and it can move in any direction, we use the relative position to encode the local layer MDP state; accordingly, the XY bounds are $[-1,1]\times[-1,1]$. As we subsequently mention in section E.2, the problem scenarios simulate routes on a grid representing $10\text{\,}\mathrm{km}$ $\times$ $10\text{\,}\mathrm{km}$, therefore velocity and acceleration limits for the agent (drone) are scaled accordingly to reflect real-world limits. All parameter values are based on the original set of values used for the DREAMR experiments, and detailed in DreamrHHP/data/paramsets of the DreamrHHP.jl repository. For example, in scale-1.toml, the XYDOT_LIM parameter which sets the speed threshold in each direction is $0.002$ representing $20\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$, which is the maximum speed of the DJI Phantom 4.

E.2 Problem Scenarios

As we mentioned in section 4.2, we use the exact same problem scenarios from the original DREAMR paper ([4], sec. V-B), and the following description is largely derived from the reference work. Since we care about the higher-level decision making framework, we abstract away physical issues like obstacles, collisions, and so on. The $1\times 1$ unit grid represents an area of $10\text{\,}\mathrm{km}$ $\times$ $10\text{\,}\mathrm{km}$ (approximately the size of north San Francisco). Each episode lasts for $30$ minutes, with $360$ timesteps or epochs of $5\text{\,}\mathrm{s}$ each. An episode starts with between $50$ to $500$ cars, with more added randomly at later epochs (up to twice the initial number). Therefore the total number of cars over the episode is $100$ to $1000$.

A new car route is generated by first choosing two endpoints more than $2\text{\,}\mathrm{km}$ ($0.2$ units) apart. We choose a number of route waypoints from a uniform distribution of $5$ to $15$, denoted $\mathrm{U}(5,15)$, and a route duration from $\mathrm{U}(100,900)$s. Each route is either a straight line or an L-shaped curve. The waypoints are placed along the route with an initial ETA based on an average car speed of up to $50\text{\,}\mathrm{m}\text{\,}{\mathrm{s}}^{-1}$. At each epoch, the car position is propagated along the route and the ETA at each remaining waypoint is perturbed with $p=0.75$ within \pm$$5\text{\,}\mathrm{s} (highly dynamic). The route geometry is simple for implementation, but we represent routes as a sequence of waypoints. This is agnostic to route geometry or topology and can accommodate structured and unstructured scenarios. In all problems, the agent begins at the centre of the grid, and the goal is near a corner.

E.3 Further Evaluation Details

We provide some additional details on the baselines, compare the various methods on two additional sets of problems of the DREAMR domain, and provide average computation times.

E.3.1 Baselines

We use the same two-level Receding Horizon Control baseline that is domain-specific to the DREAMR problem from the original paper. We used the open-source implementation from their accompanying repository (link). For UCT, table 1 lists the relevant parameters for the three variants of UCT which are used as a baseline in section 4.2. The parameters are based on the open-source implementation (link) of Double Progressive Widening that we use. We had a wider full range of parameters (14 different sets) and we chose the three most varying and representative ones. For all UCT variants, we used the default values for the double progressive widening branching parameters $(k_{a}=10,\alpha_{a}=0.5,k_{s}=10,\alpha_{s}=0.5)$

E.3.2 Two Additional Episode Sets

The two important parameters that define a DREAMR problem scenario are the number of active car routes over the episode (a measure of the size of the context set and the number of valid mode sequences to the goal) and the probability of perturbing the ETA at remaining route waypoints within \pm$$5\text{\,}\mathrm{s} (a measure of how dynamic the context set is). For the results in section 4.2, of the main paper, we generated $1000$ episodes where the number of cars is $100$ to $1000$ and the waypoint perturbation probability is $p=0.75$. This replicates the set of scenarios from the original DREAMR paper. Additionally, we generate two other sets of $1000$ episodes; the parameter values for all three sets are depicted in table 2.

For all three problem sets, the average trajectory cost for each method to reach the goal is depicted in fig. 6 and the average number of mode switches in table 3. The procedure for obtaining the statistics is the same as in section 4.2 of the main paper. For HSP we only plot the values for $\beta=0.75$ as the values for $0.55$ and $0.95$ are nearly identical to it, as was the case for Set 1. The plots demonstrate that HSP consistently outperforms the baselines .They also show how the relative performance gap between HSP and the baslines increases with more transit vehicle routes (from Set 1 to Sets 2-3), i.e. more valid mode sequences and a greater benefit in energy saved for choosing good connections and making them in time.

Furthermore, between Sets 2 and 3 (which have the same transit vehicle routes for each episode but Set 2 has a higher ETA perturbation probability), there is a slight decrease in cost incurred for most algorithms in Set 3 compared to Set 2 (as well as a slight increase in the number of mode switches). This is not unexpected as the context is changing more dynamically in Set 2, potentially increasing the number of missed connections or time spent hovering due to a delay. For HSP, the relative performance change between Sets 2 and 3 is minimal, far lower than the relative change for the other algorithms. Specifically, between Sets 2 and 3, the relative decrease in cost for the UCT variants is roughly $2.5\%$ and for RHC it is roughly $6.6\%$ (deterministic replanning is the most sensitive to perturbation), but for HSP it is only $0.2\%$. Therefore, HSP is the most robust to the variation in the context set, manifested here as the waypoint time-stamp perturbation.

E.3.3 Computation Times

For the results in section 4.2, in the main body, we focused on solution quality as our primary metric. However, we were also motivated to mitigate the increase in complexity due to considering uncertainty in our hybrid planning framework. Table 4 compares the computation time for our framework against the baseline methods. For HSP and RHC, we show planning times for both the global layer and the local layer. In practice, as we mentioned in section 3.3, after the first plan, the global layer could plan periodically in the background while the local layer executes the current plan. Therefore, in general, the decision frequency of the HSP framework is that of the local layer. Only when there is an asynchronous or event-driven interrupt, would the HSP framework be bottlenecked by the global layer. For UCT, we show the computation time required to select the next action (there are no layers).

The computation times were obtained by randomly choosing a subset of the $1000$ episodes of Set 1 (where the number of car routes is between $100$ and $1000$) and computing the average elapsed time for the various methods. The computation time depends on the exact size of the context set for the episode, so we provide an approximate range of values. Also, for the global layer, we only use the computation time for the first $25\%$ searches (subsequent global searches become trivially fast as the agent gets closer to the goal). As table 4 demonstrates, compared to RHC, which repeatedly solves a deterministic problem, we are only slightly less efficient computationally, at both global and local layers. Compared to UCT, the local layer of HSP is at least two times faster, which we do expect as the local layer is looking up a policy computed offline while UCT is doing online planning. Even the global planning of HSP, which has to search over mode sequences and transition points, is only up to two times slower than UCT.