Shortest path with acceleration constraints: complexity and approximation algorithms

We introduce a variant of the Shortest Path Problem (SPP), in which we impose additional constraints on the acceleration over the arcs, and call it Bounded Acceleration SPP (BASP). This variant is inspired by an industrial application: a vehicle needs to travel from its current position to a target one in minimum-time, following pre-defined geometric paths connecting positions within a facility, while satisfying some speed and acceleration constraints depending on the vehicle position along the currently traveled path. We characterize the complexity of BASP, proving its NP-hardness. We also show that, under additional hypotheses on problem data, the problem admits a pseudo-polynomial time-complexity algorithm. Moreover, we present an approximation algorithm with polynomial time-complexity with respect to the data of the original problem and the inverse of the approximation factor ϵ\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\epsilon$$\end{document}. Finally, we present some computational experiments to evaluate the performance of the proposed approximation algorithm.


Introduction
The Shortest Path Problem (SPP in what follows) is one of the best known within the field of combinatorial optimization. Let G = (V, A) be a directed graph and c ij be the cost of an arc (i, j) ∈ A . Let n = |V| and m = |A| . Let o, d ∈ V , o ≠ d , be an origin and a destination node, respectively, and let P od be the set of all directed paths in G from o to d. Each P ∈ P od is a subset of A made up of adjacent arcs, the first one starting at o and the last one ending at d. In the SPP the minimum cost path between o and d is searched for, that is, formally:

3
Shortest path with acceleration constraints: complexity… traverse arc (i, j) when departing from i at time t. In these problems the first-in firstout (FIFO) property is usually assumed. This states that along any arc (i, j) an earlier arrival at node j can never be attained by a later departure at node i (that is, for t ′ > t it always holds that t � + c ij (t � ) > t + c ij (t) ). Many variants of this problem are presented, for instance, in [6]. In [18], it is shown that a variant where there is a constraint on the waiting time at each node is NP-hard. Further variants are discussed in [16]. In these variants a penalty or a limit is imposed on the total waiting time spent at a given subset of the nodes. It is proved that some variants are polynomially solvable, while some others are NP-hard depending on the subset of nodes, on the fact that a penalization or a limit is imposed, and on the magnitude of the penalty parameter or of the waiting limit parameter.
Another interesting variant is SPP with time windows. In this case we associate to each node i ∈ V an interval [a i , b i ] , and to each arc (i, j) a time t ij to traverse it. Only paths in P od where each node of the path is visited within the allowed time window are feasible. The case where only elementary paths are feasible has been proven to be strongly NP-hard in [10]. The problem is still NP-hard if we remove such restriction and, thus, if we allow multiple visits to the same node. However, in such case pseudo-polynomial time algorithms have been proposed (see, for instance, [8]). A variant with additional costs associated to the departure times at the different nodes of the path is studied in [17]. In fact, a polynomial time algorithm exists if the FIFO property holds (see [9]). Moreover, if waiting is allowed at a node, then each instance for which the FIFO property does not hold can be transformed into an equivalent one for which the property holds, thus making the problem solvable in polynomial time (see [19]).
In previous work [1], we introduced a further variant of SPP called BASP (Bounded Acceleration SP) and proposed a solution algorithm. The interest for BASP arises from an industrial application, namely the optimization of automated guided vehicles (AGVs) motions in automated warehouses (see also the description in Sect. 2). Typically, AGVs follow pre-assigned routes associated to a graph, in which nodes represent operating positions and arcs represent connecting paths. AGVs motions must satisfy constraints on maximum speed, and tangential and transversal accelerations. BASP consists in finding the path in P od and the speed profile that allows travelling in minimum-time. The speed profile must satisfy maximum speed, acceleration and deceleration constraints, associated to each arc. In [1], we presented a solution algorithm (adaptive A * ) for k-BASP, a subclass of BASP. Roughly speaking, a BASP instance is a k-BASP, with k ∈ ℕ , if the maximum number of nodes of a path that can be traveled with a speed profile of maximum acceleration, followed by one of maximum deceleration, starting and ending with null speed, without violating the maximum speed constraint, is smaller than k (see [1] for details). In [1], we proved that k-BASP has polynomial time-complexity with respect to the graph size. The algorithm introduced in [1] computes the optimal trajectory between a pair of nodes and adaptively determines the value of k.
In some applications, AGVs can freely move within the facility, without having to follow any predetermined circuit. In such cases, the motion planning problem is usually addressed employing environmental representations such as cell decomposition methods. Among these, [5] presents an algorithm based on a modification of Dijkstra's algorithm, in which arcs weights depend on previously visited edges. Note that k-BASP in [1] has some similarities with the problem discussed in [5]. Indeed, in both problems, the incremental for adding an edge to a path does not depend on the complete path, but only on the k last visited nodes. However, the problem addressed by Cowlagi and Tsiotras [5] is different from BASP. Indeed, the goal of the problem in [5] is to obtain a feasible path taking into account the vehicle maximum curvature radius. On the other hand, in our work, we focus on selecting the optimal path among a set of possible paths which are already known to be feasible, while, at the same time, obtaining the optimal speed profile. Moreover, we do not assume that the incremental cost depends only on the last k visited nodes, but present conditions on problem data under which this property holds.
In this paper we discuss BASP in further detail, presenting two novel complexity results. Namely, we show that BASP is NP-hard and that, under additional assumptions, BASP admits a pseudo-polynomial time algorithm. We also introduce an -approximation algorithm based on the discretization of the admissible speeds at the nodes of the graph and show computational experiments comparing the -approximation algorithm with the solution algorithm presented in [1].
The paper is structured as follows. In Sect. 2 we describe and motivate BASP. In Sect. 3 we discuss how to compute optimal speed profiles along a single arc (with given initial and final speeds), and along a fixed path belonging to P od . In Sect. 4 we present the aforementioned complexity results. More precisely, in Sect. 4.1 we prove that BASP is NP-hard, while in Sect. 4.2 we prove that, under the assumption of integer data, BASP admits a pseudo-polynomial time algorithm. In Sect. 5 we present the -approximation algorithm. Finally, in Sect. 6 we present some computational experiments.

Problem description and motivation
In this section, we describe a new variant of SPP. The following values are associated to each arc (i, j) ∈ A: • the length ij ; • the maximum speed v max ij which can be reached along the arc; • the minimum acceleration (or maximum deceleration) a min ij < 0 and maximum acceleration a max ij > 0 along the arc.
Moreover, a zero initial speed is assigned to node o and a zero final speed is assigned to node d. We would like to select a path P ∈ P od minimizing the time needed to run along the path by fulfilling the maximum speed, the maximum and minimum acceleration constraints along the arcs, and the boundary zero conditions on o and d. Note that SPP is a special case of this problem. Indeed, SPP turns out to be equivalent to the case a max ij = +∞ and a min ij = −∞ for all arcs (i, j) ∈ A . In this case, the speed can be changed instantaneously, so that the running time along an arc is minimized by traversing it at the maximum allowed speed along the arc. Therefore, the minimum 1 3 Shortest path with acceleration constraints: complexity… time to traverse arc (i, j) is c ij = ij v max ij and the problem becomes a standard SPP with such costs c ij associated to the arcs. Since SPP corresponds to the case of unbounded acceleration limits, the proposed variant of SPP is also called Bounded Acceleration SP (BASP in what follows). In the search for an optimal solution of BASP, we should not only search for an optimal path in P od , but we should also define the speed profile along such path. As we will see later on in Sect. 3, the minimum-time speed profile along an arc (i, j) is fully determined by the initial speed v i at node i and by the final speed v j at node j. Thus, the speed v i at each node i traversed by a path P od is also part of the decision process. Such speeds must fulfill a continuity constraint, that is, if arc (i, j) is followed by arc (j, k) along the path, the final speed along arc (i, j) must be equal to the initial speed along arc (j, k). Note that the speeds at the origin node o and at the destination node d are fixed in advance.
Before proceeding, we further motivate the interest for the BASP variant of SPP. As previously mentioned, the interest comes from an industrial application. In automated warehouses, an AGV is required to pick some good up at some point of the warehouse and deliver it at some other point. The AGV moves along predefined paths and is allowed to choose among different routes at some exchange points. Formally, the exchange points as well as the points where goods are picked up and delivered represent the nodes of the graph, while the predefined routes correspond to the arcs of the graph (see, for instance, Figs. 8,13). Speed and acceleration limits differ across the different routes. For instance, if a route is a straight line, its maximum speed is higher than the maximum speed allowed along a curved route. Moreover, different speed limits may also be imposed at different points of the warehouse. For instance, if a route lies in a part of the warehouse where also human operators are working, then, for safety reasons, a lower speed limit along this route is imposed with respect to another route lying in a part of the warehouse where human operators are not allowed to work. Due to various reasons, also acceleration bounds may differ from arc to arc. For instance, in the same warehouse, flooring materials can vary from location to location. To avoid wheel slipping, we need to set maximum acceleration bounds depending on the floor frictional force, that varies according to the flooring material. Moreover, the presence of ramps along some arcs implies different acceleration and deceleration bounds. Finally, to reduce lateral oscillations, we should impose lower acceleration bounds along high-curvature connecting paths.

Remark 2.1
In the problem description we imposed a constant speed limit v max ij along each arc (i, j). In a more realistic setting the speed limit should be a function of the position along a given route: is the maximum allowed speed at position s along the route represented by arc (i, j). In particular, along a curve the speed limit varies with the curvature ray at each position along the curve. For ease of exposition, we only discuss the case of a constant speed limit along an arc (although different from arc to arc). However, at the cost of some additional technicalities, the following discussion can be extended also to the more realistic case of variable speed limits along each arc.

Remark 2.2
In a real industrial scenario, the AGV may encounter moving obstacles, such as human operators and/or other AGVs. For safety, an AGV typically halts or slows down if it perceives the presence of a human operator. When the obstacle is no longer sensed, the AGV can compute a new motion by solving a new instance of BASP, starting from its current location.
The presence of various AGVs leads to a variant of the Multi-Agent Path Finding (MAPF) problem (see, for instance, [22]), with speed and acceleration constraints. Due to the simultaneous planning of multiple AGVs, this problem is quite different from BASP (and in general much more complex), and is outside the scope of the present work. However, note that some solution approaches for standard MAPF (that does not consider velocity or acceleration bounds), such as conflict-based search, make use of sub-procedures that involve the solution of a number of standard SPP problems. Similarly, one could guess that BASP solutions could be used as basic building block to solve MAPF with speed and acceleration bounds.
In order to better appreciate the difference between SPP and BASP we can also consider the example illustrated in Fig. 1  Shortest path with acceleration constraints: complexity…

Minimum traveling time along arcs and paths
As previously pointed out, for a given arc (i, j), we are able to compute the maximum speed at which we can traverse the arc and, consequently, the minimum time needed to traverse it, as soon as we know the initial and final speeds v i and v j . Indeed, the maximum squared speed at each position s ∈ [0, ij ] along the arc is given by the following piecewise-linear function: where here and in what follows, w denotes the squared speed (so w i = v 2 i , w j = v 2 j , and so on). The result is illustrated in Fig. 2a. Starting at node i with squared speed w i , the speed is increased with the maximum possible acceleration a max ij , until the maximum allowed squared speed w max ij along the arc is reached. Such maximum speed is maintained as long as possible (null acceleration) and, finally, the speed is decreased with the maximum deceleration a min ij in order to reach squared speed w j at node j. Note that it might be possible that the maximum speed along the arc is not reachable. In such case, first the speed is increased with maximum acceleration a max ij and then decreased with maximum deceleration a min ij to reach the final squared speed w j (no constant speed portion is present in the maximum speed profile). Consequently, the minimum time to traverse arc (i, j) is the following function of the two (squared) speeds w i and w j whose solution can be derived in closed form. It is worthwhile to remark that w i and w j cannot be arbitrary values. Indeed, besides the obvious constraints w i , w j ≤ w max ij , it must also hold that If, for instance, the first inequality is not fulfilled, then it is not possible to reach the squared speed w j at node j by starting with squared speed w i at node i even by accelerating as much as possible (acceleration a max ij ). This is illustrated in Fig. 2b. In this case we set c ij (w i , w j ) = +∞.
While we reported above the formula for the minimum time needed to traverse an arc (i, j), given the initial and final squared speeds w i and w j , we further point out that a slightly more complicated formula can be derived to compute the minimum time to traverse a full path in P od . In particular, Then, it can be proved (see, for instance, [4] for a proof under more general assumptions) that the optimal (squared) speed profile is so that the minimum time to travel along the given path is The result is illustrated in Fig. 3a-c. In Fig. 3a we notice that F 0 (function F over the interval [s 0 , s 1 ] ) is obtained by starting at node o with zero speed and increasing the speed with maximum acceleration a max ij until the maximum squared speed w max ij is reached, and then keeping this speed until the end of the arc (as it is the case in the figure) or, alternatively, until the end of the arc is reached (in which case the maximum squared speed w max ij is not reached). Next, F 1 (function F over the interval [s 1 , s 2 ] ) is obtained similarly but with initial speed F(s 1 ) , while F 2 (function F over the interval [s 2 , s 3 ] ) is constant and equal to the maximum squared speed w max kh since F(s 2 ) > w max kh . In Fig. 3b we notice that B 3 (function B over the interval [s 2 , s 3 ] ) is Construction of a maximum (squared) speed profile along a path obtained in a way completely similar to F 0 . Moving backward, we start from the final node d ≡ h with squared speed w h = 0 , and we increase the speed with the maximum deceleration value a min kh until either we reach the maximum allowed squared speed w max kh along the arc, in which case we keep such speed until the beginning of the arc at node k, or we reach the beginning of the arc without reaching the maximum speed along the arc (in the figure we are in the first situation). Function B 2 (function B over the interval [s 1 , s 2 ] ) is obtained similarly but with initial speed at node k equal to B(s 2 ) . Finally, B 1 (function B over the interval [s 0 , s 1 ] ) is constant and equal to the maximum squared speed w max The optimal (squared) speed profile is illustrated in Fig. 3c and is obtained as the point-wise minimum of the functions F and B. While Figs. 2a, b and 3a, c give an intuitive illustration of the optimal speed profiles both for the case of a single arc (i, j) and for the case of a full path from o to d, we point out that these results can be derived as special cases of a more general result presented in [4] for problems with upper speed limits depending on the position along the arcs, as discussed in Remark 2.1.
We also make the following remark which states two properties of the optimal speed profile and will be useful later on.

Remark 3.1
Let W(s) be the optimal squared speed profile along a given path According to the previous discussion, the optimal speed profile along a path P from o to d with |P| + 1 nodes (here and in what follows |P| denotes the length of path P) is identified once the speeds at nodes of the path are known. Indeed, along any edge (i, j) with given initial and final squared speeds w i and w j , the optimal speed profile function is equal to (1) and the traveling time is equal to (2). Then, if we denote by P(i) ∈ V the node at position i along P, the fastest time T(P) for traversing P is the solution of the following problem:

3
Shortest path with acceleration constraints: complexity… We denote by * (P) the optimal solution of this problem and in Appendix B we will describe a recursive procedure to find it, similar to the one employed to define the optimal speed profile function W. According to the discussion above, the problem of finding the speed law which guarantees to traverse a fixed path from o to d at a minimum time, under speed and acceleration constraints, is easily solvable even in closed form. But our aim is to compute the minimum time to move from o to d by searching within all paths in P od . This is the BASP problem: As we will see in Section 4, this problem turns out to be NP-hard.

Complexity results
In this section we provide two complexity results for BASP. The first result proves NP-hardness of BASP, while the second proves that BASP admits a pseudo-polynomial time algorithm.

NP-hardness
In this section we prove that, differently from SPP, the BASP variant is NP-hard. We show this by a polynomial reduction of the NP-complete Partition problem to BASP. In the Partition problem, given a set N = {1, … , n} of positive integer values 1 , … , n , we would like to establish whether N can be partitioned into two subsets N 1 and N 2 such that Given an instance of the Partition problem, we polynomially reduce it to an instance of BASP as follows. Let G = (V, A) be such that: We set the following lengths for the arcs: For what concerns the maximum speed values, we set √ W , while we set the maximum acceleration a max ij = 1 and the minimum acceleration a min ij = −1 for all arcs except (n + 1, n + 2) , while we set a max n+1,n+2 = 0 and a min n+1,n+2 = − 1 2W . Note that, according to the imposed restrictions, a max n+1,n+2 should be strictly larger than 0. However, the result proved with null maximum acceleration can be extended, by continuity, to any sufficiently small and positive maximum acceleration. The origin node o is node 0, with zero speed, while the destination node d is n + 2 , with zero speed. An example of BASP instance derived from the Partition problem with n = 3 is illustrated in Fig. 4.
We prove the following. The proof of Proposition 4.1 is presented in Appendix A.

Remark 4.1
The complexity result given above shows that BASP is NP-hard even in case all arcs except one share the same acceleration and deceleration bounds. It is still an open question whether NP-hardness still holds if all arcs have the same bound. However, in Sect. 5.1 we will show that optimal solutions for this case are elementary paths (paths with no node repetition). This allows to derive sharper approximation results.

Pseudo-polynomial algorithm
Although BASP is NP-complete, the following proposition shows it admits a pseudo-polynomial algorithm under the assumption of integer data.
, a max is equal to 1 and a min is equal to -1 along all arcs except (4, 5) where a max 45 = 0 and a min Shortest path with acceleration constraints: complexity… Proof First, we observe that at optimal solutions there is a finite number of speeds which can be reached at each node and the squares of such speeds are integer values. Indeed, the squared speed at some node i is: • either equal to w max ij , for some j such that (i, j) ∈ A , which is an integer value by assumption; • or equal to w max ki , for some k such that (k, i) ∈ A , which is again an integer value by assumption; and squared speed at node i which is an integer value due to integrality of all the data (see Fig. 5). Note that j 1 may be the starting node o, in which case and squared speed at node i which is, again, an integer value due to integrality of all the data (see Fig. 6). Note that node j k−1 may be the destination node d, in which case w j k−1 = 0 ( j k is not included in this case). Thus, the set W of different possible squared speeds can be taken equal to the set of all integers between 0 and W = max (i,j)∈A w max ij . Now we create a new graph with node set V × W , that is, each node is a pair made up by a node in V and one of the possible squared speeds in W . Thus, the number of nodes is W |V| . For what concerns the arc set, in this graph an arc between node (i, w i ) and node (j, w j ) exists if there exists an arc (i, j) ∈ A and, moreover, c ij (w i , w j ) < +∞ , that is, there exists a feasible speed profile along arc (i, j) with initial squared speed w i and final squared speed w j . Then, the number of arcs is limited from above by W 2 |A| . The distance associated to this arc is the minimum time for a path from i to j with the boundary conditions w i and w j , which can be easily computed through (2), as discussed in Sect. 3 (recall that, in case w i and w j are not feasible, as illustrated in Fig. 2b, then the arc is removed). Then, we can solve our problem by applying, for instance, Dijkstra algorithm to this graph. Dijkstra's complexity is O(m + n log(n)) and is, thus, polynomial with respect to the size and the data of the original problem, which proves pseudo-polynomiality. ◻ Remark 4.2 While Proposition 4.2 has been proved under the assumption of integer data, it can also be extended to rational data. In such case the squared speeds which can be reached by optimal solutions are not integer values but are multiple of a rational number 1 t , where t depends on the problem data. Of course, the size of the extended graph increases with t. The approximation algorithm discussed in the following Sect. 5 is motivated by the need to consider a discretization step larger than 1 t in order to have a graph of manageable size. Shortest path with acceleration constraints: complexity…

Approximation algorithm
In this section we present an approximation algorithm for BASP with a complexity that is polynomial with respect to the size and the data of the original problem and the inverse of the approximation factor. The idea is to discretize the squared speeds in order to obtain a finite set of possible squared speeds at each node of the graph. In this way, imposing that the initial and final squared speeds along each arc belong to the discretized set of squared speeds, the set of possible speed profiles over each arc becomes finite. Hence, we can define an extended graph that enables us to solve this discretized version of the problem by means of Dijkstra's algorithm. Differently from Proposition 4.2, here arc lengths, accelerations and speed bounds need not be integer values (actually, the approach could also be extended to the case of non-constant speed bounds along arcs as discussed in Remark 2.1). In this case, we just discretize the squared speeds and impose the additional constraint that the squared speeds at the beginning and at the end of each arc belong to the set of discretized squared speeds. Let with W = max (i,j)∈A w max ij , be the set of discretized squared speeds with discretization step h. Then, The discretized problem is defined over a graph that extends graph G of the original problem. Namely, the extended graph G � = (V � , A � ) is defined as follows: where we notice that we bound from above the possible squared speeds at node i by the maximum squared speeds along the incoming and outgoing arcs of node i, while where we recall that c ij ( i , j ) = +∞ means that no feasible profile is able to travel from i to j with initial and final squared speeds equal to i and j , respectively, while c ij ( i , j ) < +∞ , as defined in (2), is the optimal travel time along arc (i, j) with the given initial and final squared speeds.

Remark 5.1 The cost for constructing the extended graph
the number of arcs of the extended graph |A ′ | is bounded from above by |A| ⋅ |Ω h | 2 and the cost for checking whether an arc exists or not in the extended graph, that is the cost for checking conditions (3), is constant.
Once the extended graph has been defined, the proposed approximation algorithm is nothing but the application of Dijkstra's algorithm to solve the discretized problem. More precisely, we search for the shortest path connecting nodes (o, 0), (d, 0) ∈ V � over graph G ′ , where the cost of arc ((i, i ), (j, j )) ∈ A � is equal to the value c ij ( i , j ) < +∞ defined in (2). Then, we have the following complexity result for the approximation algorithm.

Proposition 5.1 The complexity of the approximation algorithm is
Proof The approximation algorithm is Dijkstra's algorithm applied on the extended graph G ′ , so that its complexity is O(|A � | + |V � | log |V � |) . Then, the result immediately follows by observing that As a next step, we want to obtain an estimate of the absolute error in terms of travel time of the discretized solution returned by the approximation algorithm with respect to the continuous solution of the original BASP problem. To this end, let us consider the optimal path P * ∈ P od of the original BASP problem, with the corresponding squared speeds at each node Our aim is to build a feasible solution of the discretized problem traversing the same arcs as path P * and whose speed profile is not above the optimal speed profile of the original BASP problem but is as close as possible to it. Building such solution requires some care. It is tempting to proceed as follows: for each node i in the optimal path P * , with squared speed w i in the optimal speed profile of BASP, replace w i with that is, with the largest discretized speed bounding from below w i . Unfortunately, this does not work. Indeed, let us consider some arc (i, j) ∈ P * with the related optimal squared speeds w i and w j , and let i and j be chosen as in (7). Unfortunately, in the extended graph, arc ((i, i ), (j, j )) may not exist, since a situation like the one displayed in Fig. 2b may occur. Formally, according to (3), it may happen that either ij ij < i . Therefore, we need to proceed differently to build a solution of the discretized problem starting from the optimal solution of BASP. For x, h ∈ ℝ , h > 0 , we denote by ⟨x⟩ = max{i ∈ ℤ | ih ≤ x} the maximum multiple of h lower than or equal to x. First, we reformulate the discretized problem as follows. Let P be an assigned path from o to d. The minimum time T h (P) to traverse P in the discretized problem is: Problem (8) is obtained by adding to Problem (5) the last constraint, imposing that the coordinates of have to be multiples of h. We call h (P) ∈ ℝ |P|+1 a vector that corresponds to the solution of (8). Note that Problem (8)  Recall that this problem can be solved by the application of Dijkstra's algorithm to the extended graph described above. The following lemma compares the optimal values and solutions of Problems (5) and (8).

vector inequalities are intended component-wise).
Proof Statement (i) is a consequence of the fact that Problem (8) is obtained by adding a constraint to Problem (5). For Statement (ii), first note that, if 1 , 2 ∈ ℝ |P|+1 are feasible values for in Problem (5), then also their component-wise maximum = max{ 1 , 2 } is feasible (see, for instance, [4]). By contradiction, if h (P) ≰ * (P) then = max{ h (P), * (P)} is feasible for Problem (5). Note that the objective function (5) is strictly decreasing, that is, if 1 ≥ 2 and 1 ≠ 2 , then f ( 1 ) < f ( 2 ) . Since ≥ * (P) and ≠ * (P) , it follows that f ( ) < f ( * (P)) , contradicting the optimality of * (P) . ◻ Next, the following lemma, whose proof is given in Appendix B, gives an upper bound on the components of the difference vector * (P) − h (P) . Note that such components are nonnegative in view of part ii) of Lemma 5.1.

Lemma 5.2
The following holds for all i ∈ {1, … , |P| + 1}: In order to use Lemma 5.2, we need to find an upper bound on |P * | , the number of arcs of the optimal path. This can be done as follows.

Lemma 5.3 An upper bound for the number of arcs
|P * | of the optimal path P * is given by where • P SPP is the shortest path connecting o to d over the original graph, when the cost of each arc (i, j) ∈ A is equal to ij ; Proof A feasible solution for BASP is obtained by the arcs of path P SPP , along which, starting from null speed at node o, we first accelerate with acceleration a min until we reach speed w , then we keep the speed constant and, finally, we decelerate with deceleration a min until we reach the final null speed at node d. The time t SPP to travel along path P SPP with the given speed profile, is an upper bound for the travel time of the optimal path P * . We have that Next we need a lower bound for the time needed to travel along any arc (i, j) ∈ A . We denote this time with t min and a lower bound for it is min ∕

√W
. Then, the ratio of t SPP to t min is an upper bound for the number of arcs |P * | in the optimal path, that is, ◻ In the following, we estimate the difference T h (P h ) − T(P * ) ≥ 0 between the optimal values of the discretized BASP (9) and BASP (6) (nonnegativity follows form part i) of Lemma 5.1). Path P h , corresponding to the solution of the discretized BASP, can be different from P * and T h (P h ) ≤ T h (P * ) . Hence, T h (P h ) − T(P * ) ≤ T h (P * ) − T(P * ) . Quantities T(P * ) and T h (P * ) correspond to the optimal values of Problems (5) and (8) on the same path P * . In order to bound the difference T h (P * ) − T(P * ) , we need a further lemma, giving, for some arc (i, j), an upper bound for the difference c ij (w, z) − c ij (ŵ,ẑ) (that is, the time difference

3
Shortest path with acceleration constraints: complexity… between the times needed to travel arc (i, j) with boundary speeds w, z and ŵ,ẑ , respectively, with w ≤ŵ and z ≤ẑ ). The lemma will be proved in Appendix C.
Then, we can provide an estimate on T h (P * ) − T(P * ) (which also bounds T h (P h ) − T(P * ) ) by summing up the contributions of all arcs of P * .

Proposition 5.2 The following bound holds:
where P SPP is the shortest path connecting node o to node d.
Proof First observe that, for any arc e = (P * (i), P * (i + 1)) , it holds that as a consequence of Lemma 5.4 (first inequality) and of Lemma 5.2 (second inequality). Then, it follows that where the first inequality derives from (11), whilst the second one follows from Lemma 5.3. ◻ We can also provide an estimate of the relative error.

Proposition 5.3
Given ∈ (0, 1) , the relative error of the approximated problem with h = C , where C is a constant that depends on the problem data, is 1 + , that is, Proof Let t min be the travel time of the arc of shortest length min assuming that the squared speed along it is constantly equal to W (this is a lower bound for the travel time along any arc). Then, by Proposition 5.2, we have the following estimate over the relative error: with ◻ The following theorem states a time-complexity and an error estimate result for the approximated problem. ∈ (0, 1) , let h = C with C defined as in (12). Then,

Theorem 5.1 Given
T(P * ) ≤ 1 + , that is the solution returned by the approximation algorithm has relative error at most , and the approximation algorithm has time-complexity Proof The thesis directly follows from Propositions 5.1 and 5.3. ◻

The case of uniform acceleration bounds
As previously mentioned in Remark 4.1, it is still unclear whether the case where all arcs share the same acceleration and deceleration bounds, denoted with a max Shortest path with acceleration constraints: complexity… and a min , respectively, is NP-hard or not. However, we can prove the following result.

Proposition 5.4 If that is, all arcs have the same acceleration and deceleration bounds, then the optimal solution of BASP is an elementary path (a path that does not contain loops).
Proof Let us consider two distinct paths and ′ from node o to node d such that is elementary and ⊂ ′ . Since ′ contains the entire path , it must contain a cycle C, as shown in Fig. 7. We will show that the time needed to traverse is lower than the time needed to traverse ′ , so that we can restrict the attention to elementary paths. Let j be the first node of cycle C in ′ as indicated in Fig. 7. Let v j be the initial speed at node j along path . Along cycle C we can speed up. If, after traveling cycle C, the new speed at node j is v C > v j , then there is a time gain t in traveling the final sub-path j → d , since the new initial speed at j is higher. However, there is also a loss of time t to travel along cycle C, We prove that the latter is always greater than the former. Indeed, a lower bound t min for the time needed to travel along cycle C is obtained by assuming that the maximum acceleration a max can be kept along the cycle without hitting maximum speed bounds: Now, let be the length of the sub-path j → d . An upper bound Δt max on the time difference between the time to travel along the sub-path j → d with starting speed v j and the time to travel along the same sub-path with starting speed v C is obtained by assuming that we can keep the maximum acceleration along the sub-path: It follows that which holds true in view of v C > v j .
The same reasoning applies to maximum deceleration a min . Indeed, suppose that along path the agent has to decelerate before j to reach the speed v j . In this case, along cycle C we could speed down, allowing to maintain a higher speed along the sub-path o → j . However, it is possible to prove that the lost time t to travel along cycle C is always larger than the time gain t in traveling the initial sub-path o → j with a higher final speed, and so is faster than ′ in any case. ◻ This result is important because it allows us to state that in this sub-case the upper bound for the number of arcs of the optimal path P * depends only on the number of nodes n. In particular Therefore, we can use this bound rather than the one derived in Lemma 5.3 in all subsequent results. So, for instance, the upper bound in Proposition 5.2 can be set equal to 4n 2 h a min √w.

Computational experiments
In this section we test the approximation algorithm introduced in Sect. 5. We consider two real-life industrial scenarios of warehouses. The problem data have been provided by packaging company Ocme S.r.l., based in Parma, Italy. The simulations have been performed on a 2.7 GHz Intel Core i5 dual-core with 8GB of RAM. The acceleration and deceleration bounds on both scenarios are given by + = 0.28 m s −2 and − = −0.18 m s −2 , for arcs with mean curvature smaller than or equal to 0.25. Whilst on arcs whose mean curvature is larger than this value, the bounds are given by + = 0.14 m s −2 and − = −0.09 m s −2 . The first scenario is modeled as a graph G = (V, A) with 368 nodes and 679 edges and is depicted in Fig. 8.
The speed bounds on both scenarios are constant for each arc but vary from arc to arc, according to the associated paths curvatures.  In what follows, we compare the approximation algorithm of Sect. 5 with the adaptive A * algorithm for k-BASP presented in [1]. Note that the adaptive algorithm solves a k-BASP instance in polynomial time complexity with respect to the number of edges and vertices of the associated graph, but its complexity is exponential with respect to k. In the following, the adaptive A * algorithm for k-BASP is initialized with k = 3 . Again, we refer the reader to Ardizzoni et al. [1] for a more in depth discussion.
Given the graphs associated to the considered automated warehouses, we generated the extended graphs associated to them for different values of the discretization step h of the squared speeds. For the first scenario, step h takes values in a set H of thirty logarithmically spaced values between 0.005 and 0.5 m 2 s −2 , hence, we considered thirty different sets of discretized squared speeds Ω h and extended . We generated 1000 random pairs of source-target nodes {(s i , t i )} i∈{1,…,1000} in V × V . Then, for each of the previous pairs of source-target nodes, we considered the corresponding pair (s h h and ran Dijkstra's algorithm on G ′ h in order to obtain a trajectory starting at node s i with null speed and ending at node t i with null speed, for h ∈ H . Figure 9 shows the box-and-whisker plot of the computational times of Dijkstra's algorithm for different values of the discretization step h for solving the 1000 randomly generated instances and compares them with those of the adaptive A * algorithm for k-BASP on the same set of instances.
Note that, as the discretization step h increases, the number of discretized squared speeds decreases, hence, the number of nodes and edges in graph G ′ h decreases as well, making Dijkstra's algorithm explore a smaller graph and run faster. Also, observe that, for values of h greater than 0.015159 m 2 s −2 , the mean computational times of Dijkstra's algorithm on extended graphs G ′ h (represented by a green line with circles) are better than that of the adaptive A * algorithm for k-BASP (represented by   Figure 10 represents the box-and-whisker plot of the relative error. Note that we set a tolerance on the relative error of the trajectories obtained with the approximation algorithm, relative errors smaller than 10 −4 are not considered. This roughly corresponds to an absolute error of the order of 10 −2 s.
A good compromise for this scenario could be h = 0.063448 m 2 s −2 which is associated to a mean computational time that is 5.76 times faster than that of the adaptive A * algorithm for k-BASP, while at the same time maintaining a mean relative error of the order of 2 ⋅ 10 −2 . We could also exploit the approximation algorithm just for obtaining a path and then compute the optimal speed profile along such a path by the procedure described in Sect. 3. In this way we could employ a bigger discretization step h for achieving small computational times while maintaining high precision. The speed planning algorithm described in Sect. 3 has linear-time computational complexity with respect to the number of nodes of the path. Figure 11 shows the box-and-whisker plot of the computational times of the approximation algorithm, as in Fig. 9, summed with those of the speed planning algorithm applied on the obtained paths for the 1000 randomly generated instances on the first scenario. Figure 12 shows the box-and-whisker plot of the relative error on the travel time of the trajectories obtained coupling the approximation algorithm of Sect. 5 with the speed planning one. Again, we set a tolerance on relative errors of 10 −4 . In this case the mean relative errors are on average two orders of magnitude smaller than those presented in Fig. 10 and the percentage of solutions with a relative error smaller than 10 −4 ranges from 93.1% with h = 0.5 m 2 s −2 to 100% for h = 0.005 m 2 s −2 .  For this scenario, if we fix h = 0.5 m 2 s −2 , we get a mean computational time that is 46.35 times faster than that of the adaptive A * algorithm for k-BASP, while obtaining a solution with a mean relative error of 4 × 10 −3 , which is one order of magnitude smaller than that of the approximation algorithm alone with h = 0.063448 m 2 s −2 .
It is worthwhile to remark that there always exists a sufficiently small value h such that the optimal path of the discretized problem coincides with the optimal path of the continuous problem. Indeed, if h is chosen in such a way that the absolute error of the approximation algorithm, bounded from above as discussed in Proposition 5.2, is lower than the difference between the traveling times of the best and second-best path, then the approximation algorithm returns the best path (that is, according to the notation introduced in Sect. 5, P h is equal to P * ). However, the difference is not known in advance and even in case it were known, the choice of h based on the upper bound stated in Proposition 5.2 may lead to an excessively small value.
The second scenario is modeled as a graph with 2419 nodes and 4255 edges and is depicted in Fig. 13.  logarithmically spaced values between 0.03 and 1 m 2 s −2 , hence, we considered ten different sets of discretized squared speeds Ω h and extended graphs . As for the previous scenario, we generated 1000 random pairs of source-target nodes {(s i , t i )} i∈{1,…,1000} in V × V and tested the approximation algorithm on such pairs. Figure 14 shows the box-and-whisker plot of the computational times of Dijkstra's algorithm for different values of the discretization step h for solving the 1000 randomly generated instances and compares them with those of the adaptive A * algorithm for k-BASP on the same set of instances.
Observe that, for all values of h, the mean computational times of Dijkstra's algorithm on extended graphs G ′ h (represented by a green line with circles) are better than that of the adaptive A * algorithm for k-BASP (represented by a dashed line). However, note that the median computational time of the latter is almost the same as that of the approximation algorithm with h = 0.03 m 2 s −2 (both represented as horizontal red lines within their corresponding blue boxes). This is due to the fact that the adaptive A * algorithm for k-BASP presents a small group of outliers with very high computational times compared to its median. On the other hand, as the discretization step h increases, so does the relative error on the travel time, which, for values of h ≥ 0.45876 m 2 s −2 , is larger than 10 −1 . Figure 15 represents the box-and-whisker plot of the relative error for which we set a tolerance of 10 −4 .
A good compromise for this scenario could be h = 0.21046 m 2 s −2 which is associated to a mean computational time that is 107.3 times faster than that of the adaptive A * algorithm for k-BASP, while at the same time maintaining a mean relative error of the order of 3 ⋅ 10 −2 . Figure 16 shows the box-and-whisker plot of the computational times of the approximation algorithm, as in Fig. 14, summed with those of the speed planning algorithm described in Sect. 3 applied on the obtained paths for the 1000 randomly generated instances on the first scenario.  Figure 17 shows the box-and-whisker plot of the relative error on the travel time of the trajectories obtained coupling the approximation algorithm of Sect. 5 with the speed planning one. Again, we set a tolerance on relative errors of 10 −4 . In this case the mean relative errors are on average two orders of magnitude smaller than those presented in Fig. 15.
For this scenario, if we fix h = 0.21046 m 2 s −2 , we get a mean computational time that is 99.24 times faster than that of the adaptive A * algorithm for k-BASP, while obtaining the exact solution up to a tolerance of 10 −4 on the relative error in 97.6% of the cases.
Observe that the construction of the extended graphs can be a time-consuming operation. One could alternatively run a dynamic programming approach by generating arcs only when (and if) needed. However, the construction has to be performed only once and the extended graphs can be stored for future use. The memory occupancy of the extended graphs of the scenarios we considered varies from 36 KB for G ′ h with h = 0.5 m 2 s −2 to 258 MB for G ′ h with h = 0.005 m 2 s −2 for the first scenario, and from 78 KB for G ′ h with h = 1 m 2 s −2 to 56.1 MB for G ′ h with h = 0.03 m 2 s −2 for the second scenario.

Conclusions
Motivated by an industrial application, in this paper we addressed a variant of the Shortest Path Problem (SPP). The variant is called BASP (Bounded Acceleration SP) since speed and acceleration constraints are imposed over the arcs. Differently from SPP, where the traveling time of an arc is constant, in BASP the traveling time depends on the initial and final speed along the arc and, thus, due to the speed and acceleration constraints, it also depends on the arcs preceding and following it along a path. We proved that BASP is NP-hard but also that, under the assumption of integer data, it admits a pseudo-polynomial time algorithm. We also proposed an approximation algorithm based on the solution of an SPP problem over an extended graph. The extended graph is defined by discretizing the admissible speeds at the nodes of the graph. Finally, we performed different computational experiments on two real-life industrial scenarios in order to evaluate the performance of the approximation algorithm.

A Proof of Proposition 4.1
Each path P from node 0 to node n + 2 has the following structure with i 1 < i 2 < ⋯ < i r . Let us denote by N P = {i 1 , i 2 , … , i r } the set of intermediate nodes in P. The length of path P is W 2 + (P) , where (P) = ∑ i∈N P i is the length of the path up to node n + 1.
Before proceeding with the proof we give the intuition behind it. We will show that for each path with length (P) > W 2 up to node n + 1 , the traveling time from node 0 to node n + 2 is larger than the traveling time of a path with length (P) = W 2 up to node n + 1 (if any). This will simply follow from the fact that the former path is longer and the speed along the final arc (n + 1, n + 2) is the same in both cases. Moreover, we will show that also for a path with length (P) < W 2 up to node n + 1 , the traveling time from node 0 to node n + 2 is larger than the traveling time of a path with length (P) = W 2 up to node n + 1 (again, if any). In this case this will be proved by observing that the speed reached along this path at node n + 1 will be lower than √ W (the speed reached by any path with length (P) ≥ W 2 ). Since it is not possible to accelerate along arc (n + 1, n + 2) , the traveling time along this arc will be higher with respect to paths with length (P) ≥ W 2 up to node n + 1 . In particular, we will show that the increase of the time needed to travel along arc (n + 1, n + 2) will be larger than the time saved along the sub-path from node 0 to node n + 1 with respect to a path for which (P) = W 2 . In conclusion, we will show that the Partition problem has a yes answer if and only if the optimal value of the BASP instance is the traveling time of a path with length (P) = W 2 up to node n + 1 . In what follows we prove the result in a formal way. Let us first assume that (P) < W 2 . In this case, according to the discussion in Sect. 3, the maximum speed which can be reached at node n + 1 is v u = a max t P , where a max = 1 and t P fulfills (P) = 1 2 a max t 2 , that is, t P = Along the final arc of the path (n + 1, n + 2) the maximum acceleration is null, while the maximum deceleration is − 1 2W . Then, the optimal speed profile along this arc is obtained by keeping the speed v u = √ 2 (P) for a portion of the arc with length

3
Shortest path with acceleration constraints: complexity… W 2 − 2 (P)W , while in the last portion, with length 2 (P)W , the speed is decreased with the maximum possible deceleration − 1 2W . Now we denote by t 1 L 1 ( (P)) the time to run along the first portion of the arc, while we denote by t L 2 ( (P)) the time to run along the second part of the arc. We have that while t L 2 ( (P)) fulfills the following condition that is, Thus, the overall time to traverse such paths is Now, let us consider paths P such that (P) ≥ W 2 . A lower bound for the time needed to run along the path up to node n + 1 is given again by the solution of the following simple equation (P) = 1 2 a max t 2 , which is t P = √ 2 (P) . Note that this is a lower bound since with the maximum acceleration we would reach the speed v u = a max t P = √ 2 (P) ≥ √ W , so that we stop accelerating as soon as we reach the maximum speed √ W . Since (P) ≥ W 2 , we have that the lower bound can be further bounded from below by √ W . Finally, we observe that such lower bound can be attained if and only if (P) = W 2 , that is, if and only if the Partition problem admits a solution. Now, over the last arc (n + 1, n + 2) we can decrease the speed to 0 at node n + 2 by keeping the maximum deceleration − 1 2W over the whole arc. Then, the time needed to traverse this arc, denoted by t L 3 , fulfills so that t L 3 = 2W 3∕2 . Then, a lower bound for the time needed to traverse such paths is T 2 = √ W + 2W 3∕2 and, as already pointed out, such lower bound is attained if and only if (P) = W 2 . Figure 19 illustrates the optimal speed profiles for three distinct paths P 1 , P 2 and P 3 , fulfilling (P 1 ) < W 2 , (P 2 ) = W 2 , and (P 3 ) > W 2 , respectively. The three profiles are depicted, respectively with a dashed, a continuous and a dotted line (up to distance (P 1 ) they are overlapping). We notice that for the shortest path P 1 we are unable to reach the maximum squared speed W, so that along arc (n + 1, n + 2) we first increase the speed as much as possible with acceleration a max n+1,n+2 (in fact, we keep the speed constant since a max n+1,n+2 = 0 ), and then we decrease the speed as fast as possible with deceleration a min n+1,n+2 . For the intermediate path P 2 we reach the maximum squared speed exactly at the end of the path and then, along arc (n + 1, n + 2) we decrease the speed as fast as possible with deceleration a min n+1,n+2 . Finally, for the longest path P 3 we reach the maximum squared speed before the end of the path, and we keep the speed constant in the last part of the path, while along arc (n + 1, n + 2) we decrease the speed as fast as possible with deceleration a min n+1,n+2 . Now, if we are able to prove that T 2 < T 1 ( (P)) when (P) ≤ W 2 − 1 , then we are done. We have that By the change of variable x = √ 2 (P) , this can be rewritten as which can be easily seen to be negative for all x ≤ W 3∕2 W+1 . Now, recalling the definition of x, this means that T 2 − T 1 ( (P)) is negative if Then, the result follows by observing that

3
Shortest path with acceleration constraints: complexity… reported there, the solution of Problem (5) can be found as follows. Define vectors , , ∈ ℝ |P|+1 such that: Note that the components of vector are equivalent to the values of function W(s) defined in Sect. 3 evaluated at nodes of the path, while along each arc of the path function W(s) is defined according to (1). As a consequence of Theorem 2 of [3] the solution of Problem (13) is * (P) = . On the other hand, Problem (8) can be rewritten in the following form: Note that, strictly speaking, this is a relaxation of Problem (8), since we have not included the last constraint, namely ⟨w i ⟩ = w i , i ∈ {1, … , �P� + 1} . In fact, as shown in [3], at least one constraint is active at each component w h i (P) , i ∈ {1, … , |P| + 1} , of the optimal solution of (14). Since 0 = ⟨0⟩ and, for all x ∈ ℝ , ⟨⟨x⟩⟩ = ⟨x⟩ , necessarily ⟨w h i (P)⟩ = w h i (P) , i ∈ {1, … , |P| + 1} (which means that the optimal solution of (14) fulfills the last constraint in (8)). Also Problem (14) belongs to the class of Problem 5 in [3]. Hence, its solution can be computed with the same procedure used for Problem (13). Namely, define vectors h , h , h ∈ ℝ |P|+1 such that: Again, by Theorem 2 in [3], h (P) = h . Note that, by their definitions, h ≤ and h ≤ . Now we are ready to prove Lemma 5.2. We first prove that for all i ∈ {1, … , |P| + 1}: T h (P) = min

3
Shortest path with acceleration constraints: complexity… Remark C.1 Given c > 0 and x ∈ [0, c] , it holds that Next lemma establishes a bound from above for the error along the first part of the arc.

Fig. 18
Comparison on a generic arc between the optimal speed profile starting from ŵ and ending at ẑ and the one starting from w and ending at z Fig. 19 Maximum (squared) speed profiles along three paths P 1 , P 2 , P 3 fulfilling (P 1 ) < W 2 , (P 2 ) = W 2 , and (P 3 ) > W Lemma C.1 Given arc (i, j) and 0 ≤ w,ŵ, z,ẑ , with w ≤ŵ and z ≤ẑ the following bound holds: Proof If ŵ = 0 , then also w = 0 , hence Otherwise, if ŵ > 0 , we can bound the error over [0,s] as follows where, in the first inequality, we used the fact that, for |w −ŵ| ≤ŵ + 2a max ijs , in view of Remark C.1 it holds that: and that, for |w −ŵ| ≤ŵ,

◻
The error over the second part [s,s] of arc (i, j) can be bounded from above as follows.
Lemma C.2 Given arc (i, j) and 0 ≤ w,ŵ, z,ẑ , with w ≤ŵ , z ≤ẑ , the following bound holds:  Finally, we establish a bound for the error along the third part of the arc.

Lemma C.3
Given arc (i, j) and 0 ≤ w,ŵ, z,ẑ , with w ≤ŵ , z ≤ẑ , the following bound holds: Proof The proof of Lemma C.3 is analogous to that of Lemma C.1. ◻ Now, we can prove Lemma 5.4, that is, we can provide an estimate on the absolute error over the entire arc (i, j).