1 Introduction

Mobile robots demand more and more autonomy. They are necessary for increasingly complex tasks with limited or even nonexistent human supervision. Examples of these tasks in unstructured irregular terrains include searching for victims in emergency response scenarios [1], harvesting and/or fumigating in precision agriculture [2] or even planetary exploration [3] among others. A commonly used strategy to make mobile robots autonomously navigate is the use of path planning algorithms. These algorithms serve to calculate a path that connects the robot location to any other reachable destination [4]. Moreover, minimizing energy consumption as much as possible is desirable to make the robot perform a larger number of tasks. Thus, the path provided to the robot should be the one requiring the least amount of energy. In other words, the path planning algorithm should be optimal. Yet, this is not trivial as the idea of optimality differs between the different existing path planning approaches. We mention here three categories that are used for outdoors navigation in static environments: Sampling-based methods, Graph Search methods and Partial Derivative Equation (PDE) Solving methods [4]. These methods either create a graph that models the environment or use a pre-existing one. This graph is processed to retrieve the path from it.

Sampling-based methods model and process the graph at the same time. They iteratively create samples over a region until the initial and the goal locations are connected. This is the case for the Rapidly-exploring Random Tree (RRT), which starts from one location and creates samples (branches) resembling the growth of a tree [5]. It is asymptotically optimal: after connecting with the goal location, more samples can be created before retrieving the path. This serves to incrementally find better solutions as time runs. Nevertheless, this entails a high demand for memory resources, since usually large numbers of iterations and samples are needed to come relatively close to the global optimal solution [6].

Other approaches process a graph with pre-existing connected samples, named nodes. Graph Search methods calculate a path that is made up of these nodes. The most basic of them is Dijkstra [7]. Heuristic methods like A* evolved from Dijkstra to speed up computation [8]. Their main drawback is that the shape of the path is restricted to connecting neighbouring graph nodes. Any-angle algorithms reduce this restriction by introducing the possibility of making the path be made up of non-neighbouring graph nodes. One of them is Field-D*, an algorithm implemented on the NASA rovers for the Mars Exploration Rovers (MER) and Mars Science Laboratory (MSL) missions [9, 10]. However, the final solution is not globally optimal [11]. Algorithms like Theta* are still not globally optimal but provide better results in distance-minimizing path planning [12].

There is another group of path planning algorithms able to generate smooth paths in a globally optimal way: the Partial Derivative Equation (PDE) solving methods. They use also a pre-existing graph, but, unlike Graph Search methods, they do not compute any parent-child relationships between nodes. This removes the necessity to restrict the path to pass through intermediate nodes or edges. Instead, PDE Solving algorithms assign a characteristic direction to each node. This direction corresponds to the optimal heading of any path that crosses the node in question. This direction comes from solving the PDE in the respective node [13]. In this way, the path comes via interpolation with these characteristic directions. The equation used models an optimization formula, and its form determines which solver methods are suitable. For example, the Fast Marching Method (FMM) models the propagation of a wave over the graph [14]. The rate of propagation varies with cost values defined over each node, which can be non-uniform [15]. Nevertheless, as will be seen later, some forms of cost require the consideration of the vehicle heading. This means the cost may change with direction, i.e. it may be anisotropic. FMM produces sub-optimal results when this cost is anisotropic [13]. Another PDE Solving method, the Ordered Upwind Method (OUM), addresses this anisotropy in exchange for increasing the computational load [16]. It can hence find the optimal path considering the variations of cost according to the heading of the robot.

A path is optimal given the locomotion capabilities and the terrain. Therefore, the planner must address the terrain-vehicle interaction. The traverse of uneven surfaces is inherently anisotropic: the direction of the robot affects its energy consumption. For instance, a robot does not invest the same energy ascending and descending the same slope. Several works in the past proposed the use of Graph Search planners along with anisotropic power consumption models based on friction and gravity [17,18,19,20,21]. It is worth mentioning that while not ensuring finding the globally optimal path, these Graph Search approaches allow the use of heavy discontinuities such as forbidding ascending certain inclinations or substituting them by discontinuous zig-zag manoeuvres [20]. Other work uses nonlinear programming techniques to optimize energy while considering dynamic constraints, connecting grid nodes with feasible non-holonomic trajectories [22]. It still relies on Graph Search based methods for global planning, but this has also proven useful for local motion planning, and trajectory tracking. There is recent research oriented toward finding paths for the optimal ascension of slopes using Sampling-based methods like RRT* [23], where the difference in the slip between going straightly or diagonally through them is considered [24]. Another approach proposes the use of Sampling-based methods like SBMPO to find energy-minimizing paths considering any rise in elevation [25].

This paper presents a novel continuous and differentiable anisotropic cost model compatible with the OUM. This cost model serves to make the planner find energy-minimizing paths by addressing the energy consumption of a robot on inclined terrains. The structure of this paper comes in the following form. Section 2 presents the optimization problem that formulates the anisotropic path planner. It also presents its application on irregular terrains, focusing on the algorithm functioning and the cost function requirements. Later on, Sect. 3 fully describes the proposed anisotropic cost function. It starts with the mathematical background and later explains how this function addresses gravity, friction and slip. Next, Sect. 4 presents the results of a simulation experiment, which serves to validate its potential. Last, Sect. 5 provides our conclusions regarding the usage of an anisotropic cost function along with the OUM for navigating through scenarios containing slopes. Moreover, this section provides as well an overview of possible related future work.

2 Optimal anisotropic path planning

Fig. 1
figure 1

Different paths \(\Gamma = \Gamma ({\tilde{x}}_\textrm{o}, {\tilde{x}}_\textrm{g}, \cdot )\) connecting \({\tilde{x}}_\textrm{o}\) and \({\tilde{x}}_\textrm{g}\) can be defined upon different tangent directions \(\vec {u}(\Gamma ({\tilde{x}}_\textrm{o}, {\tilde{x}}_\textrm{g}, s))\), where s is the path length. The role of the anisotropic path planner is to find the optimal path among them

Full size image

An anisotropic PDE serves as the cornerstone to formulate and solve the optimal path planning problem. Its solution leads to finding the optimal path on an irregular surface for a mobile robot. The region \(\Omega \subset {\mathbb {R}}^2\) encloses the region of interest for the problem, depicted as a rectangle in Fig. 1. The target of the path planner is to find the optimal path connecting the goal \({\tilde{x}}_\textrm{g} \in \Omega \) and the origin \({\tilde{x}}_\textrm{o} \in \Omega \). This path is a continuous curve that can be parametrized by the path length s. In this way, the function \(\Gamma ({\tilde{x}}_\textrm{o}, {\tilde{x}}_\textrm{g}, s)\) returns the position at such a curve given a value of s, while \(\Gamma ({\tilde{x}}_\textrm{o}, {\tilde{x}}_\textrm{g}, \cdot )\) refers to the whole path. This path is the optimal one between \({\tilde{x}}_\textrm{o}\) and \({\tilde{x}}_\textrm{g}\) by complying with Eq. (1) and has \(s_\textrm{g}\) as its total length. The function \(Q(\Gamma ({\tilde{x}}_\textrm{o}, {\tilde{x}}_\textrm{g}, s), \vec {u}(\Gamma ({\tilde{x}}_\textrm{o}, {\tilde{x}}_\textrm{g}, s)))\) is the cost function. This function returns the value of cost given a location (in this case, \(\Gamma ({\tilde{x}}_\textrm{o}, {\tilde{x}}_\textrm{g}, s)\)) and a direction (in this case, \(\vec {u}(\Gamma ({\tilde{x}}_\textrm{o}, {\tilde{x}}_\textrm{g}, s))\), which is explained later). \(T({\tilde{x}}_\textrm{o}, {\tilde{x}}_\textrm{g})\) is the amount of cost accumulated along the optimal path. This amount is referred to as the total cost. It is the objective function that the planner has to minimize when searching for the optimal path.

$$\begin{aligned} \begin{aligned} T({\tilde{x}}_\textrm{o}, {\tilde{x}}_\textrm{g})&= \min _{\Gamma ({\tilde{x}}_\textrm{o}, {\tilde{x}}_\textrm{g}, \cdot ) \in \Omega } \left\{ \int _{0}^{s_\textrm{g}} Q(\Gamma ({\tilde{x}}_\textrm{o}, {\tilde{x}}_\textrm{g}, s), \right. \\&\quad \vec {u}(\Gamma ({\tilde{x}}_\textrm{o}, {\tilde{x}}_\textrm{g}, s))) \text {d}s \bigg \} \end{aligned} \end{aligned}$$
(1)

The formulation of the total cost is based on the Dynamic Programming Principle (DPP). Under this principle, the total cost between two points such as the origin \({\tilde{x}}_\textrm{o}\) and the goal \({\tilde{x}}_\textrm{g}\) is the minimum possible. This makes path planning methods working with DPP globally optimal. The DPP is expressed in Eq. (2) and explained as follows. For any point \({\tilde{x}}_{ij} \in \Omega \), the total cost from \({\tilde{x}}_\textrm{o}\) to that point, \(T({\tilde{x}}_\textrm{o}, {\tilde{x}}_{ij})\), and from that point to \({\tilde{x}}_\textrm{g}\), \(T({\tilde{x}}_{ij},{\tilde{x}}_\textrm{g})\), is equal to the total cost between \({\tilde{x}}_\textrm{o}\) and \({\tilde{x}}_\textrm{g}\), only if this point \({\tilde{x}}_{ij}\) is placed at the optimal path connecting them, \(\Gamma ({\tilde{x}}_\textrm{o},{\tilde{x}}_\textrm{g}, \cdot )\).

$$\begin{aligned} \begin{aligned} T({\tilde{x}}_\textrm{o}, x_{ij}) + T(x_{ij},{\tilde{x}}_\textrm{g}) = T({\tilde{x}}_\textrm{o},{\tilde{x}}_\textrm{g}) , \\ \forall x_{ij} \in \Gamma ({\tilde{x}}_\textrm{o},{\tilde{x}}_\textrm{g}, s) \in \Omega \end{aligned} \end{aligned}$$
(2)
Fig. 2
figure 2

Schematic of the functioning of bi-OUM. Two expanding waves start from the goal and the origin nodes (\({\tilde{x}}_\textrm{g}\) and \({\tilde{x}}_\textrm{o}\)) respectively. The loop controlling each wave is coloured in purple and orange. The red arrows indicate the characteristic direction \(\vec {u}({\tilde{x}}_{ij})\) that passes through each node \({\tilde{x}}_{ij}\)

Full size image

With regards the direction \(\vec {u}({\tilde{x}}_{ij})\), it is the direction tangent to the optimal path that passes through \({\tilde{x}}_{ij}\). This is expressed in (3). This direction is also called the characteristic direction of a location \({\tilde{x}}_{ij} \in \Omega \). In this way, given how the anisotropic cost function is defined, all the optimal paths passing through the location in question will have this same characteristic direction on it.

$$\begin{aligned} \vec {u}({\tilde{x}}_{ij}) = \frac{\text {d} \Gamma ({\tilde{x}}_\textrm{o}, {\tilde{x}}_\textrm{g}, s)}{\text {d}s}\Big \vert _{{\tilde{x}}_{ij}} \end{aligned}$$
(3)

To solve (1), it is reformulated as the Hamilton-Jacobi-Bellman (HJB) equation, following the reasoning of previous works [13, 16, 26]. This equation, shown in (4), establishes the correspondence between the anisotropic cost \(Q({\tilde{x}}_{ij}, \vec {\psi })\) and the spatial gradient of the total cost \(\nabla T\). Here, the anisotropic cost function is defined not only according to the location of any grid \({\tilde{x}}_{ij}\) but also to the heading function \(\vec {\psi }\). The latter function is the direction of the robot in the XY-plane, i.e. its heading, returned as a 2D unit vector. The direction \(\vec {\psi }\) that complies with (4) in a grid node \({\tilde{x}}_{ij}\) corresponds hence to the characteristic direction \(\vec {u}({\tilde{x}}_{ij})\) passing through it. Moreover, following the DPP in (2) and the definition of the characteristic direction in (3), the total cost in (1) is expressed from \({\tilde{x}}_\textrm{o}\) to \({\tilde{x}}_{ij}\) as well as from \({\tilde{x}}_{ij}\) to \({\tilde{x}}_\textrm{g}\).

$$\begin{aligned}{} & {} \min _{\vec {\psi }} \{ \nabla T \cdot \vec {\psi } + Q({\tilde{x}}_{ij}, \vec {\psi }) \} = 0, \ \forall {\tilde{x}}_{ij} \in {\tilde{\Omega }} \nonumber \\{} & {} \nabla T = \nabla T({\tilde{x}}_\textrm{o}, {\tilde{x}}_{ij}) = \nabla T({\tilde{x}}_{ij}, {\tilde{x}}_\textrm{g}) \end{aligned}$$
(4)

The Ordered Upwind Method (OUM) uses (4) to find the optimal path on a regular grid. According to the work of Sethian and Vladimirsky [13], this algorithm has a computational complexity that depends on the anisotropy \(\Upsilon \) of the cost function: \(O(\Upsilon n_\mathrm{{nodes}} \log (n_\mathrm{{nodes}}))\). Here, \(n_\mathrm{{nodes}}\) is the number of nodes of the grid. This anisotropy \(\Upsilon \) comes as the ratio between the highest and the lowest values of cost depending on the heading. However, for the case presented in this paper, the anisotropy varies from node to node, i.e. the anisotropy here is \(\Upsilon ({\tilde{x}}_{ij})\). Therefore, the computational complexity is expected to be variable as well, bounded by the maximum existing anisotropy in the grid \({\tilde{\Omega }}\). According to Equation (5), this anisotropy \(\Upsilon ({\tilde{x}}_{ij})\) of a node \({\tilde{x}}_{ij}\) comes as the ratio between the highest possible cost at its location and the lowest according to the heading direction \(\vec {\psi }\).

$$\begin{aligned} \Upsilon ({\tilde{x}}_{ij}) = \frac{\max \nolimits _{\vec {\psi }} \{ Q({\tilde{x}}_{ij},\vec {\psi }) \} }{\min \nolimits _{\vec {\psi }} \{ Q({\tilde{x}}_{ij},\vec {\psi }) \} } \end{aligned}$$
(5)

There is an improved version of the OUM, named the bi-directional OUM. It generates the same solution in a faster way than OUM by reducing the number of visited nodes [26]. Bi-OUM calculates both \(T({\tilde{x}}_\textrm{o}, {\tilde{x}}_{ij})\) and \(T({\tilde{x}}_{ij},{\tilde{x}}_\textrm{g})\) using two propagating waves in two parallel loops: Goal wave, which starts from \({\tilde{x}}_\textrm{g}\), and Origin wave, which starts from \({\tilde{x}}_\textrm{o}\). This exploits the DPP presented in Eq. (2). Figure 2 shows this functioning with a regular grid. Each parallel loop is coloured with a different colour, either purple (Goal wave) or orange (Origin wave). As can be checked in this figure, is not necessary that the two waves propagate until reaching the starting position of the other wave. When both waves reach an intermediate linking node \({\tilde{x}}_l\) the bi-OUM process stops. This is because, as highlighted in Eq. (2), this node complies with the DPP and there is only a single path passing through it: the optimal path. This path is calculated by following the characteristic direction that is also calculated in both waves, one from \({\tilde{x}}_\textrm{l}\) to \({\tilde{x}}_\textrm{o}\) and another from \({\tilde{x}}_\textrm{l}\) to \({\tilde{x}}_\textrm{g}\). In this way, bi-OUM achieves the same optimal solution as the OUM but invests less time than the latter. This is because the bi-directional version visits fewer nodes than the original [26].

To use the HJB equation in (4) for calculating the total cost and the characteristic direction of any node, the bi-OUM employs a semi-Lagrangian discretization. It was already used by Sethian and Vladimirsky [13] when they introduced the OUM. This discretization produces Eqs. (6) and (7) as a result. They serve to calculate, in an iterative way, the value of total cost. The Origin wave uses Eq. (6) while the Goal wave uses Eq. (7). In both expressions the update process depends on two other nodes, \({\tilde{x}}_{i'j'}\) and \({\tilde{x}}_{i''j''}\), whose corresponding values of total cost are already calculated. The main difference between (6) and (7) is that the total cost from the goal node has to consider that the anisotropic cost function multiplies the input value of characteristic direction by \(-1\). This is because the propagating wave advances in the contrary direction to the heading of the robot. The iterations end when the value of \(\epsilon \in [0,1]\) is found after converging. Equation (8), also coming from the semi-lagrangian discretization, expresses how the characteristic direction is calculated, based on the value of \(\epsilon \) found.

$$\begin{aligned} \begin{aligned} T({\tilde{x}}_\textrm{o},{\tilde{x}}_{ij})&= \min \limits _{\epsilon \in [0,1]} \left\{ Q({\tilde{x}}_{ij}, \vec {u}({\tilde{x}}_{ij})) \vert \epsilon {\tilde{x}}_{i'j'} \right. \\&\quad + \left. (1 - \epsilon ) {\tilde{x}}_{i''j''} - {\tilde{x}}_{ij}\vert + \epsilon T({\tilde{x}}_\textrm{o},{\tilde{x}}_{i'j'}) \ \right. \\&\quad + \left. (1-\epsilon ) T({\tilde{x}}_\textrm{o},{\tilde{x}}_{i''j''}) \right\} \end{aligned} \end{aligned}$$
(6)
$$\begin{aligned} \begin{aligned} T({\tilde{x}}_{ij},{\tilde{x}}_\textrm{g})&= \min \limits _{\epsilon \in [0,1]} \left\{ Q({\tilde{x}}_{ij}, -\vec {u}({\tilde{x}}_{ij})) \vert \epsilon {\tilde{x}}_{i'j'}\right. \\&\quad + \left. (1 - \epsilon ) {\tilde{x}}_{i''j''} - {\tilde{x}}_{ij} \vert + \epsilon T({\tilde{x}}_{i'j'},{\tilde{x}}_\textrm{g}) \ \right. \\&\quad +\left. (1-\epsilon ) T({\tilde{x}}_{i''j''},{\tilde{x}}_\textrm{g}) \right\} \end{aligned} \end{aligned}$$
(7)
$$\begin{aligned} \begin{aligned} \vec {u}({\tilde{x}}_{ij})&= \frac{ \epsilon {\tilde{x}}_{i'j'} + (1 - \epsilon ) {\tilde{x}}_{i''j''} - {\tilde{x}}_{ij} }{|| \ \epsilon {\tilde{x}}_{i'j'} + (1 - \epsilon ) {\tilde{x}}_{i''j''} - {\tilde{x}}_{ij} \ ||} \end{aligned} \end{aligned}$$
(8)
Fig. 3
figure 3

Update process of a Considered node \({\tilde{x}}_{ij}\). Its values of total cost and characteristic direction are computed taking into account AcceptedFront nodes within the distance \(\xi ({\tilde{x}}_{ij})\) expressed in Eq. (9), also called Near-AcceptedFront nodes

Full size image

Each node has two state functions: \(S^\textrm{o}_{ij}\) and \(S^g_{ij}\). Each of these functions indicates the state of the node \({\tilde{x}}_{ij}\) according to the Origin wave or the Goal wave respectively. \(S^o_{ij}\) and \(S^\textrm{g}_{ij}\) present each one state out of the four presented in Table 1, which will be further explained later. The nodes \({\tilde{x}}_{i'j'}\) and \({\tilde{x}}_{i''j''}\) that affect the explained update processes of the total cost and the characteristic direction must present an AcceptedFront state. Besides, they must be located close enough to \({\tilde{x}}_{ij}\), with a proximity distance lower than \(\xi ({\tilde{x}}_{ij})\). This distance threshold is defined in Eq. (9). It comes in function of the anisotropy \(\Upsilon ({\tilde{x}}_{ij})\) that is present at such node, already expressed in (5). \(\xi ({\tilde{x}}_{ij})\) is also proportional to the grid resolution \(\Lambda \). Figure 3 shows how this distance works.

$$\begin{aligned} \xi ({\tilde{x}}_{ij}) = \Lambda \Upsilon ({\tilde{x}}_{ij}) \end{aligned}$$
(9)

The process followed by bi-OUM to visit the grid nodes, represented in Fig. 2, is also presented as pseudo-code in Algorithm 1. First of all, both \(S^\textrm{o}_{ij}\) and \(S^\textrm{g}_{ij}\) are initialized to Far for every node. Thereafter, each loop sets the total cost of the respective starting node to zero, a null value to the characteristic direction (since initial and final desired orientations are not defined) and the Accepted state to the respective state function. Next, each process updates the neighbours of the corresponding starting node thanks to the updateNeighbours() function. Two functions are iteratively executed by the Goal wave and the Origin wave. On the one hand, updateNeighbours() performs three actions. First, it checks whether an AcceptedFront node no longer has Far or Considered neighbours and changes its state to AcceptedInner. Second, it converts those Far nodes that are neighbours to a target node \({\tilde{x}}_{t}\) into Considered. This target node is initially the starting node from which the wave expands. Third, the total cost of the Considered neighbours to the target node is updated using either Eqs. (6) or (7), as well as the characteristic direction using Eq. (8). Thereafter, the following target node is chosen using the getNextNode() function. This function takes the existing Considered node with the lowest value of total cost, \({\tilde{x}}_{t}\), and changes its state into AcceptedFront.

Table 1 Possible states of grid nodes in the calculation of the values of total cost and characteristic direction
Full size table
figure a
Fig. 4
figure 4

Graphical representation of the slope model used to build up the anisotropic cost function. This model encompasses multiple variables such as the slope gradient \(\alpha \), the aspect \(\vec {\gamma }_{ij}\), the normal vector \(\vec {\nu }_{ij}\), the heading direction \(\vec {\psi }_{ij}\) and the relative angle \(\beta (\vec {\psi }, \vec {\gamma }_{ij})\). All of them are marked in this conceptual depiction. a 3D Perspective view of slope model. b Lateral and top views of the slope model

Full size image

Finally, the checkFinCondition() function evaluates \({\tilde{x}}_\textrm{t}\) to check if its state is AcceptedInner or AcceptedFront for both loops. If so, the waves stop and \({\tilde{x}}_\textrm{t}\) becomes \({\tilde{x}}_\textrm{l}\), indicated in Fig. 2 as a blue dot. Later on, the getPath() function returns the optimal path from two portions: one from \({\tilde{x}}_\textrm{l}\) to \({\tilde{x}}_\textrm{o}\) and another from \({\tilde{x}}_\textrm{l}\) to \({\tilde{x}}_\textrm{o}\). Each waypoint is calculated as indicated in (10), where the step distance is \(d_\textrm{step}\), and the characteristic direction for each waypoint is interpolated from the values calculated on the nodes. The final path \({\tilde{\Gamma }} = \left\{ {\tilde{\Gamma }}_\textrm{o}, {\tilde{\Gamma }}_1, {\tilde{\Gamma }}_2... {\tilde{\Gamma }}_\textrm{g}\right\} \) is obtained by joining both portions (see Fig. 2).

$$\begin{aligned} \begin{aligned} {\tilde{\Gamma }}_{k-1}&= {\tilde{\Gamma }}_{k} - d_\textrm{step} \vec {u}({\tilde{\Gamma }}_{k}) \ , \ if \ Origin \ wave \\ {\tilde{\Gamma }}_{k+1}&= {\tilde{\Gamma }}_{k} + d_\textrm{step} \vec {u}({\tilde{\Gamma }}_{k}) \ , \ if \ Goal \ wave \end{aligned} \end{aligned}$$
(10)

3 Anisotropic cost function for inclined surfaces

This section presents the steps followed to make an anisotropic cost function \(Q({\tilde{x}}_{ij},\vec {\psi })\) be based on slope parameters and compatible with the (bi-)OUM.

3.1 Slope-based anisotropy

Figure 4 shows the slope model used to represent the interaction between the robot and inclined terrain. In the absence of any kinematic configuration capable to reconfigure itself [27], the pose of the robot body will change according to the shape of the terrain surface. This change will be determined by the contact points between the robot and the surface. Nevertheless, to avoid making the formulation more complex, a simplification is made: the robot-terrain interaction is modelled after a single contact point. This simplification assumes the robot body is always parallel to an imaginary inclined plane. The normal vector of this plane, named \(\vec {\nu }_{ij}\) in Fig. 4, will be hence coincident with the Z-axis of the robot local reference frame. This imaginary plane is inclined at a certain angle from the horizontal XY plane (the plane perpendicular to the gravity vector). The value of this angle corresponds to the slope gradient \(\alpha _{ij}\) and is equal to the angle between the normal vector \(\vec {\nu }_{ij}\) and the global Z-axis as showcased in Fig. 4. The direction the slope faces, i.e. the direction in which the steepest descent occurs, is the aspect or \(\vec {\gamma }_{ij}\). It can be obtained from projecting the normal vector of the slope, \(\vec {\nu }_{ij}\), onto the 2D XY plane and normalizing it.

To make the anisotropic cost function \(Q({\tilde{x}}_{ij},\vec {\psi })\) compatible with the (bi-)OUM, it is necessary to understand its requirements. First of all, this function must always return a positive and non-zero value. Complying with this condition avoids producing undesirable local minimum points that compromise the optimality of the resulting path, making it sub-optimal [13]. Next, the solution that bi-OUM produces is viscous, meaning it not only exists but is also unique and stable. This is thanks to the fact that the computed solution omits any discontinuity present in the real solution of the total cost [16]. However, the computed solution is guaranteed to be unique if, and only if, it is ensured that the inverse of the cost function, \(1/Q({\tilde{x}}_{ij},\vec {\psi })\), is fully differentiable and convex [13, 16]. For this reason, as shown in Fig. 5 the inverse of the real cost of the robot, \(1/{\tilde{Q}}({\tilde{x}}_{ij},\vec {\psi })\), is approximated with a closed conic curve. This curve is \(1/Q({\tilde{x}}_{ij},\vec {\psi })\), which is formulated as a displaced ellipse due to its simplicity. Equation (11) expresses the polar form of this ellipse, whose radius is \(1/Q({\tilde{x}}_{ij},\vec {\psi })\). This form uses \(\beta (\vec {\psi }, \vec {\gamma }_{ij})\), which represents the angle between the robot heading \(\vec {\psi }\) and the aspect \(\vec {\gamma }_{ij}\). In this way, when \(\vec {\psi }\) coincides with \(\vec {\gamma }_{ij}\) then \(\beta (\vec {\psi }, \vec {\gamma }_{ij})\) returns a value of zero. Figure 5 depicts how the radius of this ellipse, which corresponds to the inverse of the cost \(1/Q({\tilde{x}}_{ij},\vec {\psi })\), varies with \(\beta (\vec {\psi }, \vec {\gamma }_{ij})\).

$$\begin{aligned} \begin{aligned} 0 = \frac{ [\cos _\beta \sin _\beta ]^T \left[ \begin{array}{cc} - C^\textrm{a}_{ij} C^\textrm{d}_{ij} &{} 0 \\ 0 &{} - {C^l_{ij}}^2 \end{array} \right] \left[ \begin{array}{c} \cos _\beta \\ \sin _\beta \end{array} \right] }{Q({\tilde{x}}_{ij},\vec {\psi })^2} \ + \\ + \ \frac{\left[ \begin{array}{cc} C^\textrm{a}_{ij} - C^\textrm{d}_{ij}&0 \end{array} \right] \left[ \begin{array}{c} \cos _\beta \\ \sin _\beta \end{array} \right] }{Q({\tilde{x}}_{ij},\vec {\psi })} + 1 \ , \ \ \ \beta = \beta (\vec {\psi }, \vec {\gamma }_{ij}) \end{aligned} \end{aligned}$$
(11)
Fig. 5
figure 5

Elliptical inverse of the anisotropic cost function \(Q({\tilde{x}}_{ij}, \vec {\psi })\). It serves to approximate the inverse of the real cost of the robot, \({\tilde{Q}}({\tilde{x}}_{ij}, \vec {\psi })\), while complying with the requirements of the OUM (convex, continuous, closed). The robot has a heading direction \(\vec {\psi }\) that has an angle \(\beta (\vec {\psi }, \vec {\gamma }_{ij})\) with respect to the aspect \(\vec {\gamma }_{ij}\)

Full size image

The shape of the ellipse varies with the slope gradient \(\alpha _{ij}\) through three functions: the ascent cost \(C^\textrm{a}_{ij}\), the lateral cost \(C^l_{ij}\) and the descent cost \(C^\textrm{d}_{ij}\). Each of them correspond to what the real cost \({\tilde{Q}}({\tilde{x}}_{ij},\vec {\psi })\) returns at certain values of \(\beta (\vec {\psi },\vec {\gamma }_{ij})\): \(\pm \pi \) (12), \(\pm \pi / 2\) (13) and 0 (14).

$$\begin{aligned} C^\textrm{a}_{ij}= & {} {\tilde{Q}}({\tilde{x}}_{ij},\vec {\psi } \mid \beta (\vec {\psi },\vec {\gamma }_{ij}) = \pm \pi ) \end{aligned}$$
(12)
$$\begin{aligned} C^l_{ij}= & {} {\tilde{Q}}\left( {\tilde{x}}_{ij},\vec {\psi } \mid \beta (\vec {\psi },\vec {\gamma }_{ij}) = \pm \frac{\pi }{2}\right) \end{aligned}$$
(13)
$$\begin{aligned} C^\textrm{d}_{ij}= & {} {\tilde{Q}}({\tilde{x}}_{ij},\vec {\psi }\mid \beta (\vec {\psi },\vec {\gamma }_{ij}) = 0) \end{aligned}$$
(14)

The terms \(\cos _\beta \) and \(\sin _\beta \) can be substituted by the expressions in (15), given the fact that the ellipse in Fig. 5 is symmetrical in the axis where \(\vec {\gamma }_{ij}\) is located. With these expressions in mind, the terms in (11) can be rearranged to leave this equation in the explicit form shown in (16).

$$\begin{aligned}{} & {} \begin{aligned}&\vert \cos _{\beta }\vert = \vec {\psi } \cdot \vec {\gamma }_{ij} \ \ , \ \ \vert \sin _{\beta }\vert = \vert \vert \vec {\psi } \times \vec {\gamma }_{ij} \vert \vert \end{aligned} \end{aligned}$$
(15)
$$\begin{aligned}{} & {} Q({\tilde{x}}_{ij},\vec {\psi }) =\nonumber \\{} & {} \qquad \sqrt{\left( \frac{C^\textrm{a}_{ij} + C^\textrm{d}_{ij}}{2} \right) ^2 \left( \vec {\psi } \cdot \vec {\gamma }_{ij} \right) ^2 + \left( C^l_{ij} \ \vert \vert \vec {\psi } \times \vec {\gamma }_{ij} \vert \vert \right) ^2}\nonumber \\{} & {} \qquad -\frac{C^\textrm{a}_{ij} - C^\textrm{d}_{ij}}{2} \ \vec {\psi } \cdot \vec {\gamma }_{ij} \end{aligned}$$
(16)

3.2 Energy minimizing cost function

The next step is to re-define \(Q({\tilde{x}}_{ij},\vec {\psi })\) by means of the robot energy consumption according to its dynamics. Equation (17) serves to model this consumption when climbing any slope. It is based on a model created in previous work to represent the current consumption of a robot with n wheels and wheel radius r [28]. Here we assume that the same voltage is supplied to all motors, so the current consumption is proportional to the power. The motor torque constant is referred to as \(K_m\), which serves to translate the torque into the current.

$$\begin{aligned} I(\rho _{ij}, \sigma _{ij}, \alpha _{ij}) = n K_m r \frac{m g (\rho _{ij} \cos _{\alpha _{ij}} - \sin _{\alpha _{ij}})}{\cos _{\alpha _{ij}} (1 - \sigma _{ij})} \end{aligned}$$
(17)

The terms that appear in the fraction are explained as follows. The numerator is the Drawbar Pull Resistance Force. It represents the resulting force that drags the robot. It comes as a function of the specific resistance coefficient \(\rho _{ij}\), the mass of the robot m and the gravity acceleration g. Rowe and Ross [17] use a similar expression to model the same. The coefficient \(\rho _{ij}\) estimates the ratio between the normal force, generated by the terrain surface, and the drawbar pull, produced by making the wheels roll.

Two terms are in the denominator. First, \(\cos _{\alpha _{ij}}\) arises to account for the fact that we are solving a two-dimensional path planning problem. In other words, the 2.5D elevation map is projected onto the 2D plane and the aforementioned robot speed v takes different values in the 2D projection when climbing or descending through a slope, i.e. changing its Z-coordinate. This can be understood better by checking on Fig. 4a, where the blue circle on the slope takes the form of an ellipse in its projection onto the XY-plane. Second, the slip ratio \(\sigma _{ij}\) is the difference between 1 and the ratio between \(v_{ij}\) and the estimated speed that is commanded to the wheels. In other words, \(\sigma _{ij}\) takes a value of zero when \(v_{ij}\) and the commanded speed are the same. Its value gets close to one as the commanded speed increases.

The functions \(C^\textrm{a}_{ij}({\tilde{x}}_{ij})\), \(C^l_{ij}({\tilde{x}}_{ij})\) and \(C^\textrm{d}_{ij}({\tilde{x}}_{ij})\) are defined using (17) to make the anisotropic cost function \(Q({\tilde{x}}_{ij},\vec {\psi })\) consider the path planning criterion of electric charge minimization. They are expressed in Eqs. (18), (19) and (20), respectively. All of them consider the robot speed \(v_{ij}\).

$$\begin{aligned} C^\textrm{a}_{ij}= & {} \frac{I(\rho _{ij}, \sigma _{ij}, \alpha _{ij})}{v_{ij}} \end{aligned}$$
(18)
$$\begin{aligned} C^l_{ij}= & {} \frac{I(\rho _{ij}, \sigma _{ij}, 0)}{v_{ij}} \end{aligned}$$
(19)
$$\begin{aligned} C^\textrm{d}_{ij}= & {} {\left\{ \begin{array}{ll} n K_m r \frac{m g R_\textrm{b}({\tilde{x}}_{ij})}{v_{ij} (1 - \sigma _{ij})} \ , \\ \ \alpha _{ij} \in (\alpha ^\textrm{zero}_{ij} - \alpha _{\Delta },\alpha ^\textrm{zero}_{ij} + \alpha _{\Delta }) \\ \left| \frac{I(\rho _{ij}, \sigma _{ij}, -\alpha _{ij})}{v_{ij}} \right| \ , \\ otherwise \end{array}\right. } \end{aligned}$$
(20)
Fig. 6
figure 6

Use of Bezier curve function \(R_\textrm{b}({\tilde{x}}_{ij})\) to comply with the non-zero positive cost specification in a smooth way. \(\rho _{ij} = 0.3\), \(\alpha _{\Delta } = {15}^{\circ }\) and \(\alpha ^\textrm{zero}_{ij} = \arctan _{\rho _{ij}}\)

Full size image

The function \(C^l({\tilde{x}}_{ij})\), expressed in (19), takes a value of zero for the \(\alpha _{ij}\) input. This is because, whenever the robot goes in a direction perpendicular to the aspect \(\gamma _{ij}\), the value of elevation does not change, i.e. the robot does not ascend or descend. On the other hand, a special treatment is given to \(C^\textrm{d}({\tilde{x}}_{ij})\), expressed in (20). The effect of gravity pulls the robot and reduces its energy consumption when it descends. Here, it is assumed the vehicle cannot recharge itself, so the energetic cost should always be present in the form of a positive value. Besides, it may be also desirable to prevent the robot from braking when descending through slopes, so the loss of energy in the form of heat is avoided [17]. When the robot descends, having a value of \(-\alpha _{ij}\) as input that acknowledges this robot pose, the current function from (17) could return zero or even negative values. This is incompatible with the OUM and bi-OUM requirements. Therefore, this situation is dealt with by using the Bezier function \(R_\textrm{b}({\tilde{x}}_{ij})\) that is present in (20). This function acknowledges the angle of slope gradient in which the robot would start gaining energy, \(\alpha ^\textrm{zero}_{ij}\). From this slope gradient, the robot starts braking to keep its velocity, avoiding any acceleration. \(R_\textrm{b}({\tilde{x}}_{ij})\) preserves continuity and smoothness while making the cost \(C^\textrm{d}({\tilde{x}}_{ij})\) return always positive values that increase with the slope gradient \(\alpha _{ij}\). This is showcased in Fig. 6.

Fig. 7
figure 7

Selection of region from the DEM of a crater. The elevation used is half of the original. Preparation of the map used for the simulation tests with the anisotropic planner. It is based on the shape of a real crater on the Martian surface. The resolution of the DEM is 1 m

Full size image

Given the expression in (17), the angle \(\alpha ^\textrm{zero}_{ij}\) would be equal to \(\arctan _{\rho _{ij}}\). The absolute value is used to define \(C^\textrm{d}({\tilde{x}}_{ij})\) in (20) so this function increases the cost with values of slope gradient higher than \(\alpha ^\textrm{zero}_{ij} + \alpha _{\Delta }\). The use of an absolute value was used by Rowe and Ross [17] to penalize paths that required the robot to brake and not accelerate. Here, \(\alpha _{\Delta }\) is a custom configurable angle margin. \(R_\textrm{b}({\tilde{x}}_{ij})\) marks three points using this margin: at \(\alpha ^\textrm{zero}_{ij} - \alpha _{\Delta }\), at \(\alpha ^\textrm{zero}_{ij}\) and at \(\alpha ^\textrm{zero}_{ij} + \alpha _{\Delta }\). A Bezier curve is the basis of \(R_\textrm{b}({\tilde{x}}_{ij})\) to not only preserve continuity but also smoothness while penalizing braking.

4 Experiments

This section presents an evaluation of the anisotropic cost function \(Q({\tilde{x}}_{ij},\vec {\psi })\), presented in Sect. 3 for path planning purposes. In planetary exploration missions, rovers may drive not only on horizontal surfaces but also on inclined ones. For instance, the Spirit rover was commanded over weeks to climb the Husband summit on Mars and later descend it [29]. For this reason, we prepared and carried out a numerical simulation using the model of a martian unstructured environment containing a crater. This paper focuses on executing only the bi-OUM algorithm to produce a series of paths in the mentioned scenario. This is because a past comparative study already demonstrates how bi-OUM generates better solutions and even in a faster way than others compatible with anisotropic cost functions, like OUM, RRT* and Genetic Algorithms [26].

The purpose of this experiment is twofold. First, it is of interest to analyse how the terrain affects the performance of the anisotropic cost function to find energy-minimizing paths. Second, to figure out how significant is the use of an anisotropic function in contrast with an isotropic function. The latter is usable by the Fast Marching Method (FMM), another PDE Solving method introduced in Sect. 1 that employs much lower computational complexity: \(O(n_\mathrm{{nodes}} \log (n_\mathrm{{nodes}}))\) [30]. All the code used to produce the paths of all the tests is written using the Python language and is available online.Footnote 1

The simulation test was carried out using the Digital Elevation Model (DEM) of a crater. This DEM is based on the shape of a real crater located close to where the Spirit rover landed on Mars. The data was obtained from the High Resolution Imaging Science Experiment (HiRISE)Footnote 2 repository and was adapted to make it present slopes up to 20 degrees. The main reason to do this is that the original presents pronounced slopes that would be non-traversable, and the key point in this simulation test is to have a large variety of traversable inclined surfaces. Figure 7 shows the \(80\times 80\) m portion of elevation data that was extracted from the original DEM. Figure 8 presents the elevation data describing the shape of this crater together with some contour lines to ease the visualization. The slope gradient (maximum inclination) is shown in Fig. 9 for all the points on the map.

Fig. 8
figure 8

Elevation map

Full size image
Fig. 9
figure 9

Slope gradient map

Full size image

First of all, we analyse how the specific resistance \(\rho _{ij}\) and the slip ratio \(\sigma _{ij}\) influence the energy consumption estimation provided by the anisotropic cost function. It is worth mentioning that, in a similar way to isotropic cost functions, the constant parameters that multiply the whole function, acting as gains, do not affect the location of the waypoints of the resulting path. These parameters only modify proportionally the values of total cost assigned to the nodes. Since the value returned by the specific resistance \(\rho _{ij}\) can be a real number higher than 0 and lower than 1, the simulation test is performed using a discrete set of constant values of \(\rho _{ij}\): 0.15, 0.3, 0.45, 0.6, 0.75 and 0.9. In other words, the planner is executed several times, using one out of the six different values of \(\rho _{ij}\) each time for all nodes \({\tilde{x}}_{ij} \in {\tilde{\Omega }}\).

Fig. 10
figure 10

Slip ratio function of the two models (Wheel and Track) found in the literature [31], used in the numerical simulation tests

Full size image

With regards to the slip ratio \(\sigma _{ij}\), it is defined after two models constructed from experimental data and introduced by Sutoh et al. [31]. These models were constructed from two locomotion subsystems: one using wheels and the other using tracks. These two mechanisms were simulated on an inclined sandbox, making them work with lunar regolith simulant [32]. Figure 10 depicts the two functions of slip ratio \(\sigma _{ij}\) that depend on the slope gradient \(\alpha _{ij}\). Moreover, they are expressed in Equation (21). With the increase of this slope gradient, both slip ratio functions return higher values for both models, but at different rates. The slip ratio of the Wheel model increases faster than that of the Track model.

$$\begin{aligned} \sigma _{ij} = {\left\{ \begin{array}{ll} 0.07 e^{0.1 \alpha _{ij}} &{}\quad \text {if} \ \text {Wheel} \ \text {Model} \\ 0.04 e^{0.07 \alpha _{ij}} &{}\quad \text {if} \ \text {Track} \ \text {Model} \end{array}\right. } \end{aligned}$$
(21)

In this way, there are six constant values to define the specific resistance \(\rho _{ij}\) and two functions that depend on \(\alpha _{ij}\) to define the slip ratio \(\sigma _{ij}\). Moreover, it is here created an isotropic cost function that takes the highest value of cost from the anisotropic cost function as shown in (22).

$$\begin{aligned} C({\tilde{x}}_{ij}) \equiv \max \limits _{\vec {\psi }} \{ Q({\tilde{x}}_{ij},\vec {\psi }) \} \end{aligned}$$
(22)
Fig. 11
figure 11

Values of cost returned by each anisotropic cost function. Left 3d plot of the values returned by the cost function, where two horizontal axes correspond to the ascent-descent and lateral directions, while the vertical axis is the slope gradient. Middle Polar plot of the inverse to the cost function as in Fig. 5, given several values of slope gradient. Right directional cost functions and anisotropy. The latter uses the red axes

Full size image
Fig. 12
figure 12

3d view of the scene with the resulting paths. Results from the first test using anisotropic and isotropic cost functions. The resulting paths connecting two locations, \({\tilde{x}}_\textrm{o}\) to \({\tilde{x}}_\textrm{g}\), are depicted. The origin is located at (10 10) m while the goal is at (55 50) m. Note that the 3d view is rotated to provide a better perspective of the obtained paths. The grid \({\tilde{\Omega }}\) used is a hexagonal regular one, with a resolution \(\Lambda \) of 0.5 m

Full size image

For a better understanding of how the anisotropic and isotropic cost functions are defined in this simulation test, Fig. 11 is provided. It contains 4 subfigures with 3 plots each. The left plot depicts the anisotropic cost function using a 3d representation. This representation resembles a vertical cylinder: the base axes serve to construct a polar plot while there is also a vertical axis pointing upwards (the slope gradient \(\alpha _{ij}\)). The middle plot depicts the inverse of this cost function using a 2d polar plot. In this case, the information about the slope gradient \(\alpha _{ij}\) is only represented by the use of different colours, as the vertical axis from the previous plot is projected. This view is similar to the one shown in Fig. 5. The right plot serves to represent the directional cost functions \(C^\textrm{a}_{ij}({\tilde{x}}_{ij})\), \(C^l_{ij}({\tilde{x}}_{ij})\) and \(C^\textrm{d}_{ij}({\tilde{x}}_{ij})\) as well as the anisotropy \(\Upsilon ({\tilde{x}}_{ij})\) (in red).

Figure 11a, b corresponds to the case where \(\rho = 0.15\) and the slip ratio model is the Wheel one. The difference between them is that the first one uses the anisotropic cost function while the second uses an isotropic one. As can be checked, the cost of the isotropic is equal for all directions of the robot, while in the anisotropic case there are differences according to this direction. For values of \(\alpha _{ij}\) close to 20 degrees, the inverse of the anisotropic cost takes the shape of an elongated ellipse, while the isotropic cost remains as a circle. In the isotropic case, this circle becomes smaller as the slope gradient increases, while in the anisotropic case this inverse shrinks in the ascent-descent directions, where the relative angle \(\beta (\vec {\psi }, \vec {\gamma }_{ij})\) between the heading \(\vec {\psi }\) and the aspect \(\vec {\gamma }_{ij})\) directions is either 0 or 180 degrees. The last two subfigures show the anisotropic cost function with \(\rho _{ij} = 0.3\) and each of them uses a different slip model. Figure 11c shows the use of the Wheel model, while Fig. 11d shows the use of the Track model. As can be checked, the cost in the first one is higher in the ascent and descent functions (see right images) due to the higher values of slip ratio (see Fig. 10), while the lateral cost remains the same. This also creates a different anisotropy for values of slope gradient close to 20. The difference in anisotropy is more significant by comparing the first and third subfigures, Fig. 11a, c, where an increase in \(\rho _{ij}\) entails higher anisotropy.

With the defined anisotropic and isotropic cost functions based on different configurations of terramechanic functions, a series of paths are planned. Figure 12 depicts these paths connecting two locations of interest: the origin \({\tilde{x}}_\textrm{o} = x_\textrm{o}\) and the goal \({\tilde{x}}_\textrm{g} = x_\textrm{g}\). As can be seen, those paths created with low values of \(\rho _{ij}\), between 0.15 and 0.3, go through different places than the rest. The isotropic cost functions with low \(\rho _{ij}\) get further from the slopes and surround the crater, taking more distance to reach the goal. The paths with low \(\rho _{ij}\) and generated using an anisotropic cost function traverse laterally the slopes, keeping the same elevation. This is because although the ascent and descent costs are high, the lateral cost is still lower (see Fig. 11a). To do a deeper insight into how the total cost function is calculated, Figs. 13 and 14 depict the solution calculated by two configurations: one anisotropic and one isotropic. In Fig. 13 the anisotropy makes the wave propagation from the origin enter the crater. This is because the planner acknowledges that the descent cost is cheaper than the lateral and ascent costs as also seen in the third row of Fig. 11. On the contrary, the isotropic cost from (22) prevents the wave from propagating towards the crater slopes, as it does not address the differences in cost according to direction. For this reason, in the anisotropic case (with \(\rho _{ij} = 0.3\) and using the Wheel model) the path traverses the crater, while the isotropic planner finds another path that circumvents the crater by sticking to horizontal surfaces as much as possible. Therefore, the consideration or not of the anisotropy entails different results. Figures 15 and 16 serve to highlight the implications of these differences.

Fig. 13
figure 13

Total cost calculated using anisotropic cost with \(\rho _{ij} = 0.3\) and the Wheel slip model

Full size image
Fig. 14
figure 14

Total cost calculated using isotropic cost with \(\rho _{ij} = 0.3\) and the Wheel slip model

Full size image

Figure 15 indicates the increase in the number of node updates when using anisotropic cost in contrast with using isotropic cost. This metric gives an idea about the extra computational load when calculating the solution using anisotropy. As can be denoted, the use of anisotropic cost functions entails a higher of visits to the nodes. This fact goes in line with the computational complexities of anisotropic planners like OUM and isotropic planners like FMM. As a reminder, the computational complexity of OUM is proportional to the overall anisotropy [13]. Figure 15 indicates that for low values of specific resistance \(\rho _{ij}\) the number of updates is around 20 times higher than in the isotropic case. This is because the anisotropy increases when the specific resistance is lower, especially when it gets close to zero. This can be checked by comparing the anisotropy plot of Fig. 11a (\(\rho _{ij}\) = 0.15) and Fig. 11c (\(\rho _{ij} = 0.3)\).

Fig. 15
figure 15

Increase in the number of visits to nodes when planning using anisotropic cost instead of isotropic cost. Each visit is an update of total cost and characteristic direction values using Eqs. (6), (7) and (8)

Full size image
Fig. 16
figure 16

A measure of how significant is to consider the directional dependency of the cost according to the specific resistance \(\rho _{ij}\) and to each slip ratio model. This is calculated according to Eq. (23) given the resulting anisotropic and isotropic total cost values

Full size image

Although the increase in computational load is significant, the reduction in total cost produced by considering anisotropy still needs to be measured. Figure 16 shows this reduction for all values of specific resistance and both slip models, calculated using the expression shown in (23).

$$\begin{aligned} \begin{aligned} \text {Reduction} [\%]&= \frac{T_\textrm{aniso}({\tilde{x}}_\textrm{o}, {\tilde{x}}_\textrm{g}) - T_\textrm{iso}({\tilde{x}}_\textrm{o}, {\tilde{x}}_\textrm{g})}{T_\textrm{iso}({\tilde{x}}_\textrm{o}, {\tilde{x}}_\textrm{g})} \\ {}&\quad \times 100 \end{aligned} \end{aligned}$$
(23)

The total cost of paths created using anisotropic cost is \(T_\textrm{aniso}({\tilde{x}}_\textrm{o}, {\tilde{x}}_\textrm{g})\), while \(T_\textrm{iso}({\tilde{x}}_\textrm{o}, {\tilde{x}}_\textrm{g})\) is the total cost of paths created with isotropic cost. As can be checked in Fig. 16, for \(\rho _{ij} = 0.15\) the reduction is close to \(-13\%\) for the Wheel model and \(-15\%\) for the Track model. For values of specific resistance \(\rho _{ij}\) between 0.15 and 0.45, there is some significant difference between the models. The reduction in the Track model decreases rapidly until being around \(-2.5\%\). On the contrary, the Wheel Model at \(\rho _{ij} = 0.3\) experiences a reduction of around \(-20.0\%\).

The exposed data evinces that the use of isotropic cost functions introduces undesired extra total cost in all cases. However, this extra total cost is only significant for \(\rho _{ij} < 0.45\) with the Wheel model and for \(\rho _{ij} < 0.3\) with the Track model. For this reason, the use of anisotropy in the cost function to address uneven surfaces is only recommended in those circumstances, especially when using a slip model with a high slip ratio such as the Wheel model. The user should evaluate whether the reduction in total cost compensates for the increase in the number of node updates, also quite high for lower values of specific resistance as shown in Fig. 15. For example, for the case in which this planning is carried out offline, the trade-off in question may not be a problem. Yet, the high computational load may be intractable for online planning in the onboard computer of a robotic vehicle.

5 Conclusion

In this paper, we have presented a cost function model aimed at performing anisotropic path planning on terrains containing slopes. By defining this model as the inverse polar function of a displaced ellipse, we make it compatible with the use of anisotropic PDE path planners like the bi-OUM. This cost function considers the energy consumption of the robot according to its heading when it experiences inclination. To better understand the use of the cost function, we present in this paper the results from a simulation test. These results have demonstrated in which situations the anisotropy may be beneficial for making a robot optimally traverse scenarios with slopes. In particular, two slip models were used, in which one of them, the Wheel model, was affected by the inclination more than the other, the Track model. The terrain was considered by not only its shape, through the use of the slope gradient, but also by the use of the specific resistance parameter. Since PDE planners like FMM can use isotropic functions with low computational complexity, they were used to contrast anisotropic ones. The results indicate that the higher the effect of the slip and the lower the value of specific resistance motivate the use of an anisotropic cost function instead of an isotropic one. The benefit of opting for anisotropic cost functions is more prominent in the Wheel model, given a value of specific resistance lower than 0.45.

Finally, we foresee the continuation of this work by studying the use of other PDE path planning methods such as the Fast Sweeping Method (FSM), compatible with anisotropic cost functions [33] and with turning constraints [34]. An improvement to the proposed anisotropic cost function would be the preservation of the stability, as addressed in other work [26], as well as the consideration of non-traversable areas [35]. Another would be the consideration of different functions for the slip ascending and descending as in other approaches [36]. Finally, initial and goal orientations could be addressed by adapting previous approaches done with FMM [37].