{'fraction_non_alphanumeric': 0.03554079302470417, 'fraction_numerical': 0.011604733236454225, 'mean_word_length': 4.260565687452222, 'pattern_counts': {'":': 0, '<': 7, '<?xml version=': 0, '>': 6, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A Nonlinear Constrained Optimization Framework for Comfortable and Customizable Motion Planning of Nonholonomic Mobile Robots - Part II

In this series of papers, we present a motion planning framework for planning comfortable and customizable motion of nonholonomic mobile robots such as intelligent wheelchairs and autonomous cars. In Part I, we presented the mathematical foundation of our framework, where we model motion discomfort as a weighted cost functional and define comfortable motion planning as a nonlinear constrained optimization problem of computing trajectories that minimize this discomfort given the appropriate boundary conditions and constraints. In this paper, we discretize the infinite-dimensional optimization problem using conforming finite elements. We describe shape functions to handle different kinds of boundary conditions and the choice of unknowns to obtain a sparse Hessian matrix. We also describe in detail how any trajectory computation problem can have infinitely many locally optimal solutions and our method of handling them. Additionally, since we have a nonlinear and constrained problem, computation of high quality initial guesses is crucial for efficient solution. We show how to compute them.


Introduction
In the first paper of this series (Gulati et al., 2013), we formulated motion planning for a nonholonomic mobile robot moving on a plane as a constrained nonlinear optimization problem. This formulation was motivated by the need for planning trajectories that result in comfortable motion for human users of autonomous robots such as assistive wheelchairs (Fehr et al., 2000;Simpson et al., 2008) and autonomous cars (Silberg et al., 2012). We reviewed literature from robotics and ground vehicles and identified several properties that a trajectory must have for comfort. We then reviewed motion planning literature in robotics (Latombe, 1991;Choset et al., 2005;LaValle, 2006LaValle, , 2011a and concluded that most existing methods have been developed for robots that do not transport a human user and have not explicitly addressed issues of comfort and customization. We then posed motion planning as a constrained nonlinear optimization problem where the objective is to compute trajectories that minimize a discomfort cost functional and satisfy appropriate boundary conditions and constraints. The discomfort cost functional was formulated based on comfort studies in road and railway vehicle design. We performed an in-depth analysis of conditions under which the cost-functional is mathematically meaningful, a thorough analysis of boundary conditions, and formulated the constraints necessary for motion comfort and for obstacle avoidance. We showed that we must be able to impose two kinds of boundary conditions. In the first kind, the problem is set in the Sobolev space of functions whose up to second derivatives are square-integrable. In the second kind, we must allow functions that are singular at the boundary (with a known strength) but still lie in the same Sobolev space in the interior.
The above optimization problem is infinite dimensional since it is posed on infinite dimensional function spaces. This means that we must discretize it as a finite dimensional problem before it can be solved numerically. Keeping the problem setting and requirements mentioned above in mind, it is natural to use the Finite Element Method (FEM) (Hughes, 2000) to discretize it. In this paper, we present the details of this discretization and show how to use an appropriate finite dimensional subspace for both kinds of boundary conditions. We also show the sparsity structure of the Hessian of the global problem. Some of the discretized inequality and inequality constraints, if computed naively, lead to a dense global Hessian. We avoid this and keep the global Hessian sparse by introducing auxiliary variables.
A good initial guess is crucial for solving a nonlinear optimization problem. We describe a method for computing high quality initial guesses for the optimization problem.
The choice of finite-element method for discretization, the choice of appropriate finite dimensional subspace for different types of boundary conditions, the choice of appropriate variables to be solved for, and a method to compute good initial guess together result in a fast and reliable solution method.

Background
For completeness, we briefly describe the assumptions and problem formulation that are presented in more detail in Part I of this series.

Motion of a nonholonomic mobile robot moving on a plane
We model the robot as rigid body moving on a plane subject to the following nonholonomic constrainṫ x sin θ −ẏ cos θ = 0. (1) Here dot, (˙), represents derivative with respect to t. A motion of such a body can be specified by specifying a travel time τ and a trajectory r(t) for t ∈ [0, τ ]. The orientation θ(t) can be computed from Equation (1). Essentially, θ(t) = arctan2(ṙ(t)). Ifṙ(t) is zero, which means the velocity is zero, then this equation cannot be used. If the instantaneous velocity is zero at t = t 0 , and non-zero in a neighborhood of t 0 , then θ(t 0 ) can be defined as a lim t→t0 arctan2(ṙ(t)).
Let the geometric path on which the robot moves be an arc-length parameterized curve r(s) (see Figure 1). The tangent and normal vectors to the curve are given by The signed curvature κ(s) is defined as where θ(s) is the tangent angle. From review of comfort studies in road and railway vehicles (Suzuki, 1998;Jacobson et al., 1980;Pepler et al., 1980;Förstberg, 2000;Chakroborty and Das, 2004;Iwnicki, 2006), we concluded that for motion comfort, it is necessary to have continuous and bounded acceleration along the tangential and normal directions. It is possible that the actual values of the bounds on the tangential and normal components are different. It is also desirable to keep jerk small and bounded.
We model discomfort as the following cost functional J: Here τ is the total travel time and r is the position of robot at time t ∈ [0, τ ]. ... r represents the jerk. ... r · T and ... r · N are the tangential and normal components of jerk respectively. The weights (w T and w N ) are non-negative known real numbers. We separate tangential and normal jerk to allow a choice of different weights (w T and w N ).
The weights serve two purposes. First, they act as scaling factors for dimensionally different terms. Second, they determine the relative importance of the terms and provide a way to adjust the robot's performance according to user preferences. For example, for a wheelchair, some users may not tolerate high jerk and prefer traveling slowly while others could tolerate relatively high jerks if they reach their destination quickly. We determine the the typical values of weights using dimensional analysis (see Part I (Gulati et al., 2013)).

Parameterization of the trajectory
The trajectory r can be parameterized in various ways. We have found that expressing the trajectory in terms of speed and orientation as functions of a scaled arc-length parameter leads to relatively simple expressions for all the remaining physical quantities (such as accelerations and jerks).
Let u ∈ [0, 1]. The trajectory is parameterized by u. The starting point is given by u = 0 and the ending point is given by u = 1. Let r = r(u) denote the position vector of the robot in the plane. Let v = v(u) be the speed. Both r and v are functions of u. Let λ denote the length of the trajectory. Since only the start and end positions are known, λ cannot be specified in advance. It has to be an unknown that will be found by the optimization process.
Let s ∈ [0, λ] be the arc-length parameter. We choose u to be a scaled arc-length parameter where u = s λ . The trajectory, r(t), t ∈ [0, τ ] is completely specified by the trajectory length λ, the speed v = v(u), and the orientation or the tangent angle θ = θ(u) to the curve. λ is a scalar while speed and orientation are functions of u. These are the three unknowns, two functions and one scalar, that will be determined by the optimization process.
With this parameterization, the discomfort cost functional of Equation (4) can be written as The first integral (J τ ) is the total time, the second integral (J T ) is total squared tangential jerk, and the third integral (J N ) is total squared normal jerk.
The discomfort J is now a function of the primary unknown functions v, θ, and a scalar λ, the trajectory length. All references to time t have disappeared. Once the unknowns are found via optimization, we can compute t using the following expression.
The position vector r(u) can be computed via the following integrals.
If θ(u) is known, r(u) can be computed from Equation (7). If v(u) and λ are known, t(u) can be computed from Equation (6). Using these two, we can determine the function r(t), t ∈ [0, τ ].
The expressions for tangential acceleration a T and normal acceleration a N are Here N is the direction normal to the tangent (rotated π 2 anti-clockwise). The signed curvature is given by The angular speed ω is given by

Function spaces for v and θ
We now define the function spaces to which v and θ can belong so that the discomfort J in Equation (5) is well-defined (finite). We have two distinct cases depending on whether the speed is zero at an end-point on not.
Let Ω = [0, 1] and H 2 (Ω) be the Sobolev space of functions on Ω with square-integrable derivatives of up to order 2. Let f : Ω → R. Then We can show that if v, θ ∈ H 2 (Ω), then the integrals of squared tangential and normal jerk are finite. Using the Sobolev embedding theorem (Adams and Fournier, 2003) it can be shown that if f ∈ H 2 (Ω), then f ∈ C 0 (Ω) and by extension f ∈ C 1 (Ω). Here C j (Ω) is the space of functions on Ω whose up to j th derivatives are bounded and continuous. Thus, if v, θ ∈ H 2 (Ω), then all the lower derivatives are bounded and continuous. Physically this means that quantities like the speed, acceleration, and curvature are bounded and continuous − all desirable properties for comfortable motion (Section 2.2).
We also need that the inverse of v be integrable so that J τ is finite. Inverse of v is integrable if v is uniformly positive in [0,1]. This is trivially true if v is uniformly positive, that is, v ≥ v > 0 for some constant positive v throughout the interval [0, 1]. However, v can be zero at one or both end-points because of the imposed conditions (see Section 2.5).

Positive speed on boundary
Consider the case that v is positive on both end-points. We make the justifiable assumption that the trajectory that actually minimizes discomfort will not have a halt in between. Thus, if v > 0 on end-points, it remains uniformly positive in the interior and the discomfort is finite.

Zero speed on boundary
Consider the case in which v(0) = 0. The case v(1) = 0 can be treated in a similar manner. If v(0) = 0, 1 v must not blow up faster than 1 u p where p < 1 to keep J τ finite. Thus, if zero speed boundary conditions are imposed, we will have to choose v outside H 2 (Ω). In such a case, at u = 0, it is sufficient that v approaches zero as u p where 3 5 < p < 1 or p = 1 2 . For the right end point, where u = 1, replace u with (1 − u) in the condition. v ∈ H 2 (Ω) in the interior.

Boundary conditions
The expression for the cost functional J in Equation (5) shows that the highest derivative order for v and θ is two. Thus, for the boundary value problem to be well-posed we need two boundary conditions on v and θ at each end-point − one on the function and one on the first derivative.
We now relate the mathematical requirement on v and θ boundary values above to expressions of physical quantities. We do this for the starting point only. The ending point relations are analogous.

Positive speed on boundary
First, consider the case when v > 0 on the starting point. The speed v needs to be specified, which is quite natural. The u-derivative of v, however, is not tangential acceleration. The tangential acceleration is the t-derivative and is given by Equation (8). It is vv λ .
Here v is known but λ is not. Thus specifying tangential acceleration gives us a constraint equation and not directly a value for v (0). This is imposed as an equality constraint. Similarly, fixing a value for θ on starting point is natural. We "fix" the values of θ (0) by fixing the signed curvature κ = θ λ . As before, this leads to an equality constraint relating θ (0) and λ if κ = 0. Since choosing a meaningful non-zero value of κ is difficult, it is natural to impose κ = 0. In this case θ (0) = 0 can be imposed easily. A non-convex star-shaped obstacle is shown with its "center" {x0, y0} and a distance function ρ = ρ(φ). The distance function gives a single point on the boundary for φ ∈ [0, 2π]. The robot trajectory must lie outside the obstacle.

Zero speed on boundary
If v(0) = 0, then, as seen in Section 2.4, v(u) must behave like u p for 3 5 < p < 1 or p = 1 2 near u = 0 and v (u) ∼ u q for − 2 5 < q < 0 or q = − 1 2 respectively. In this case is the appropriate strength of the singularity for v (u) at the starting point.

Obstacle avoidance constraints
We model the "forbidden" region formed by the obstacles as a union of star-shaped domains with boundaries that are closed curves with piecewise continuous second derivative. A set in R n is called a star-shaped domain if there exists at least one point x 0 in the set such that the line segment connecting x 0 and x lies in the set for all x in the set. Intuitively this means that there exists at least one point in the set from which all other points are "visible". We will refer to such a point x 0 as a center of the star-shaped domain.
Let an obstacle be specified by its boundary in polar coordinates that are centered at r 0 = {x 0 , y 0 }. Each φ ∈ [0, 2π) gives a point on the boundary using the distance ρ(φ) from the obstacle origin. The distance function ρ must be periodic with a period 2π. See Figure 2.
Suppose we want a point r = {x, y} to be outside the obstacle boundary. Define C(r) as where the subscript 2 refers to the Euclidean norm. C(r) ≥ 0 ⇐⇒ the point r is outside the obstacle. The constraint function C(r) is piecewise differentiable for all r except at a single point r = r 0 . If r = r 0 by chance, which is easily detectable, we know that the r is inside the obstacle and can perturbed to avoid this undefined behavior. Note that C(r) remains bounded inside the obstacle.

The infinite-dimensional optimization problem
We now summarize the nonlinear and constrained trajectory optimization problem taking into account all input parameters, all the boundary conditions, and all the constraints. This is the "functional" form of the problem (posed in function spaces). We will present an appropriate discretization procedure valid for all input combinations in the Section 4.
Minimize the discomfort functional J, where given the following boundary conditions for both starting point and ending point • position (r 0 , r τ ), • orientation (θ 0 , θ τ ), • signed curvature (κ 0 , κ τ ), • tangential acceleration (a T ,0 , a T ,τ ), and constraints on allowable range of • tangential acceleration (a T ,min , a T ,max ), • normal acceleration (a N ,min , a N ,max ), • angular speed (ω min , ω max ), • curvature, if necessary (κ min = 0, κ max ), and • number of obstacles N obs , and • an initial guess for (v(u), θ(u), λ), in u ∈ [0, 1], • weights w T > 0 and w N > 0. 7 The constraint on starting and ending position requires that Staying outside all obstacles requires that As a post-processing step, we compute time t as a function of u using and convert all quantities (v, θ, r, and their derivatives) from u domain to t domain.

A finite element discretization of the infinite-dimensional optimization problem
We first discuss the case when there is no singularity in the speed v. This is the case when the given boundary speeds are positive. In this case, v ∈ H 2 (Ω) as described in Section 2.4, where Ω = [0, 1]. We assume that θ ∈ H 2 (Ω) always (whether v is singular or not) because it is sufficient for the discomfort to be finite. Thus, to discretize the problem, it is natural to use the basis functions in C 1 (Ω), the space of functions that are continuous and have continuous first derivatives. This makes the second derivative of v and θ discontinuous but its square is still integrable.

Basis functions
We minimize the discomfort and satisfy all the constraints in a finite dimensional subspace of C 1 (Ω). We make the following choice.
Here χ i (u) ∈ C 1 (Ω) are basis functions for this problem. The symbol h traditionally denotes a measure of the "mesh width" to distinguish the approximate solution from the "exact" infinite dimensional solution.
The unknown scalar values α v i and α θ i for i = 1, . . . , q are the degrees of freedom (DOFs) which are to be determined via "solving" the optimization problem of Section 3. For reasonable choices of χ i (u), as q increases the finite dimensional solution approaches the exact solution.
Usually, one chooses χ i (u) ∈ Ω that are easy and cheap to compute, and have local support. Piecewise polynomial functions that have sufficient differentiability are good candidates. By having local support, we mean the functions are non-zero only over a limited interval and not on whole Ω. This has two main advantages. First, while performing integration to compute J, the product of χ i and χ j is zero over most of the interval. This leads to O(n) rather than O(n 2 ) interactions. Second, the global Hessian matrix is sparse rather than being dense. Typically for 1D problems, the matrix has a small constant band-width.
For our problem, we first divide the interval [0, 1] into n equal-size intervals of length h = 1 n . Each of these intervals is an element. Thus we have n elements and n + 1 equidistant points or nodes. The collection of elements and nodes is the finite element mesh.
At each node we define two piecewise polynomial functions that are non-zero only on the two elements surrounding the node. This is the standard cubic Hermite basis for problems posed in the H 2 space in 1D. See Figure 3 for both the functions and their first and second derivatives. The first kind of basis function is 1 on the node with which it is associated and has zero derivative there. The second function is zero on the corresponding node and has unit derivative there. Additionally, both are zero with first derivatives also zero on two surrounding nodes. Thus, in all, we have 2(n + 1) basis functions and unknown scalars each for v and θ. Because of the above-mentioned properties each basis functions belongs to C 1 (Ω) ⊂ H 2 (Ω). Note that the basis functions on the boundary points are slightly different. They are not extended outside the interval. Figure 4 shows 5 basis functions of each kind for n = 4. As seen, the boundary basis functions are truncated. It would not matter anyway what their values outside the interval are since no integration is performed outside. The other basis functions are translations of each other.

Element shape functions
In FEM practice, we define basis functions in terms of "shape functions". The shape functions are defined only on a single reference element and multiple shape functions placed on neighboring elements are joined to create a single basis function. This requires that shape functions have appropriate values and derivative on reference element boundary so that the basis functions are valid for the problem.

Singular shape functions at boundary
For the case when either one or both the boundary points have zero speed specified, v(u) is singular at the corresponding boundary with v (u) infinite with singularity u −1/3 when acceleration is also zero and u −1/2 if acceleration is non-zero. Thus, v as a function does not belong to H 2 (Ω). However, v does belong to H 2 in the interior as shown in Section 2.4.
After a FEM mesh is decided, we take this into account and do not use the above-mentioned regular shape functions for v on the element(s) near the boundary with zero speed. For the interior elements, however, no change is done and the regular shape functions are used.
We now derive the new singular shape functions for the left boundary (u ∈ [0, h]) element and the singular shape functions for the right boundary element can be derived using symmetry.
We need at least two shape functions so that the two shape functions coming from the [h, 2h] element can be matched. Denote them by ψ L 1 and ψ L 2 . The function value and the function derivative both must be matched Figure 3: Cubic Hermite basis functions and their first and second derivatives. Two kinds of basis functions are defined at each node such that each is zero everywhere except on the two elements sharing the node. The first kind, χ k,1 has a value of 1 and a zero derivative at the associated node k. Its value and derivatives are zero on both adjacent nodes. The second kind, χ k,2 has a value of zero and a unit derivative at the associated node k. Its value and derivatives are zero on both adjacent nodes. Both χ k,1 and χ k,2 are square integrable on u ∈ [0, 1].  at h. For a reference element shape function, this means ψ L 1 (1) = 1, ψ L 1 (1) = 0, ψ L 2 (1) = 0, ψ L 2 (1) = 1. Both ψ L 1 and ψ L 2 must be zero on u = 0 because the speed is zero there. Thus, the Dirichlet boundary condition is imposed explicitly. Finally, as x → 0, one of the functions must behave as x p , for p = 2 3 or p = 1 2 , to match 11 Figure 6: Singular shape functions on boundary elements for zero speed and zero acceleration boundary conditions. Two singular shape functions are defined on a boundary element. Consider zero speed condition on left element. The first shape function ψ L 1 has a value of 1 and zero derivative on the right to match the shape function φ1 of Figure 5 on the next element. The second shape function ψ L 2 has a value of 0 and a unit derivative on the right to match the shape function φ2 of Figure 5 on the next element. Both ψ L 1 and ψ L 2 have a value of zero on the left because speed is zero. The singular shape functions for zero speed boundary conditions on right are similarly defined. the singularity in Equation (13). It can be seen that the choice satisfies all these requirements. Figure 6 shows the four shape functions, two for the left element and two for the right element for the case p = 2 3 . Figure 7 shows that ψ L 1 and ψ R 1 are singular when approaching the boundary (shown for p = 2 3 only). The other two functions ψ L 2 and ψ R 2 are smooth.

The finite-dimensional optimization problem
With the choice of basis functions described above, we can express v h (u) and θ h (u) as: ψ L 1 and ψ R 1 tend to infinity as x − 2 3 as they approach the boundary. ψ L 2 and ψ R 2 are smooth.
where v i , v i , θ i , and θ i for i = 1, . . . , n are the unknown nodal values and χ i,1 (u) and χ i,2 (u) are the two kinds of basis functions described in the previous section.
For optimization, the values of cost, its gradient and Hessian, the values of constraints, and the gradient and Hessian of each constraint are required. For efficiency, it is desirable that cost and constraint Hessians be sparse. We will see later that for the Hessian of obstacle avoidance constraints to be sparse, it is useful to introduce 2N additional unknowns in the form of position Thus, our objective now is to determine the values of these unknowns and the unknown path length λ that minimize the cost functional and satisfy the boundary conditions and constraints described in (Section 3).

Numerical integration for computing the integrals in the cost functional and constraints
We use Gauss quadrature formulas to compute the integrals in the cost functional and constraints. When using m integration points in an interval, the formulas are accurate for polynomials of degree up to 2m − 1.
For the regular C 1 basis functions, which have maximum polynomial degree 3, it is easily seen that the square tangential jerk is a polynomial of degree 23, and the squared normal jerk is a polynomial of degree 17. See Equation (5). These are polynomials in u and not t. Hence, 12 Gauss points will give exact integrals up to floating point accuracy. Of course, the integrands being polynomials, the integrals corresponding to J T and J N can be evaluated without using Gauss points (if one is ready to work with complex algebraic expressions). But the other integrals, for J τ and those relating r to θ (Equation (7)), must be evaluated numerically. Hence, we use 12 Gauss points to evaluate all integrals.

Imposing constraints
We need to impose multiple equality and inequality constraints while minimizing the cost functional. Some of the equality constraints affect a single DOF each and hence they can be used to eliminate the particular unknown. The others relate multiple DOFs and must be imposed as an equality explicitly. The equality constraints are described below.
• Fix end-point positions (r 0 , r τ ) by eliminating the unknowns x and y on the first and last nodes.
• Fix end-point orientations (θ 0 , θ τ ) by eliminating the unknown θ on the first and last nodes.
• Fix end-point speeds (v 0 ≥ 0, v τ ≥ 0) by eliminating the unknown v on the first and last nodes.
• If speed on an end-point is positive, tangential acceleration a T must be specified at that end-point.
Impose vv λ = a T on that end point. Otherwise, this constraint will be automatically imposed by using singular shape functions.
• Impose specified end-point curvature κ by imposing θ = λκ on each end-point.
The inequality constraints that are not related to obstacles avoidance are as follows. We must maintain • tangential acceleration in [a T,min , a T,max ], • normal acceleration in [a N,min , a N,max ], and • angular velocity in [ω min , ω max ].
• curvature in [κ min = 0, κ max ], Note that these must be maintained for each u ∈ [0, 1] in the infinite dimensional optimization problem. For the discretized version, we choose the Gauss integration points and impose that these quantities remain in the specified range only on those points. Thus, for a mesh with n elements, each inequality above results in 12n constraints (assuming 12 points are used as discussed in Section 4.4). The values of these physical quantities on each Gauss point is a function of the DOFs on element nodes. Thus, these constraints are local. They are not affected when DOFs of non-element nodes change.
There are two important reasons to keep the values within range on Gauss points as opposed to on some other, say, uniform set of points. First, since we use v at the Gauss point to compute the integrals, it is more important that v remain non-negative there to avoid problems of large negative values of J. Second, since v and θ are already computed there it saves extra computation.

Obstacle avoidance constraints
Staying outside obstacles, if present, requires additional inequality constraints. For this we pick N uniformly separated points in the interval [0, 1] and impose the constraint that each of r on the N points remain outside each of N obs obstacles. This leads to N ×N obs constraints. In our implementation we make N = nM +(n+1), so that if the distribution is uniform, each element has M such points in the interior and each node is a point too. Two of these N points are the boundary points which must be outside all obstacles for the optimization problem to have a feasible solution. If the robot boundary is not circular and we choose P points on the robot's boundary, then the number of constraints is N × N obs × P .
We come back to obstacle related constraint relating a single obstacle and a single point on the trajectory. One could simply relate the position at the point with θ(u) and λ using Equation (7), and use Equation (14) to impose conditions on θ(u) and λ. However, because of the structure of Equation (14) and because r(u) depends on all θ DOFs of nodes that are before u, the Hessian of this constraint is not sparse. This would lead to efficiency problems when doing iterations in the numerical optimization process. Even computing the dense Hessian would be very costly as N obs and N increase. We must work around this elimination approach of imposing the obstacle related constraints.
To avoid the dense Hessian of obstacle constraint, we make two changes to the simplistic approach. First, we do not use Equation (15) for eliminating r but keep r as an unknown function. Second, we relate adjacent r's via Equation (7) as follows.
Here j goes from 1 to N − 1. What this change does is that, as long as adjacent r(u j ) and r(u j−1 ) belong to a maximum two adjacent elements, the equation above relates only a small number of local DOFs. Secondly, since each r(u j ) is now a legitimate unknown, it can be used to impose the inequality constraint Equation (14) without θ being involved.
This new approach does have a price, however. We have increased the number of unknowns and hence increased the size of the gradient vector and Hessian matrix. But this is a small price to pay considering that the sparsity is still maintained, the amount of computation does not grow, and equations are local in nature. We explore the sparsity pattern more in Section 4.7 ahead.

Element and global gradient and hessian
An important choice in the FEM discretization of any variational problem is the ordering of all the unknowns when forming the global Hessian matrix. A good choice simplifies the assembly process as well as could lead to useful structural sparsity.
We have four kinds of DOFs. For simplicity, we discuss the regular case and where boundary conditions are not yet imposed. The singular case differs in minor details only that does not affect the ordering process. The four kinds are as follows.
• four unknowns each on n + 1 node −v, v , θ, θ • N x and N y unknowns • a scalar unknown λ The unknowns are ordered in the same sequence shown above starting from u = 0 and going to u = 1.
The ordering chosen above means that each DOF except λ interacts with DOFs on two elements. The scaled arc-length parameter, λ, is global by its nature and interacts with all other DOFs. Hence, the global Hessian matrix is sparse. Some interactions lead to linear equations, so they do not affect the Hessian. This is the case for x, y, and λ interactions in Equation (20).
We now describe the structure of the finite dimensional optimization problem using a small mesh with n = 3 elements and N = 10 points for obstacle constraints as shown in Figure 8(a). In FEM, each element provides a small Hessian, typically dense, that relates all the DOFs present in that element. We have eight v and θ DOFs on each element except for the singular corner elements that have six. Figure 8(b), shows the global connectivity structure of the problem after boundary conditions are imposed on boundary v, θ, x, and y DOFs. These are marked A in Figure 8(a). The three element matrices are added to their appropriate positions. The DOFs marked B (for θ ) are constrained via equality constraints. The DOFs marked C (for v ) are constrained via equality constraints if speed is non-zero. Otherwise, it is infinite and is taken care using singular elements. The x and y DOFs do not enter the expression of cost, hence all corresponding rows and columns are empty (zero). The Hessian of obstacles constraints does contain non-zero 2 × 2 blocks relating x and y of the same point.

Convergence on decreasing mesh size
To analyze the effect of mesh size on convergence, we construct two examples of straight line motion, each with start and end positions as {0, 0} and {10, 0}, and start and end orientations as 0. In the first example, speed at both ends is 0. In the second, speed at both ends is 1. Tangential acceleration at both ends is 0. We vary the number of elements from 2 to 128 in multiples of 2 and each time solve the discomfort minimization problem for one initial guess. As the number of elements increases, the optimum cost found by the optimization process decreases. This is natural because increasing the mesh size means we're minimizing a function in a superset of degrees of freedom. As the number of elements increases, the relative change in minimum cost decreases. We compare all costs with the cost corresponding to 128 elements. We see that the 32 elements give a cost that within 0.01% of cost for 128 elements (when v > 0 on end-points). The curve for v = 0 shows a lower convergence rate and we believe that the reason behind this is using standard Gauss-Legendre quadrature for the singular elements. A more precise procedure would use specially designed quadrature scheme keeping in mind the form of singularity at end-points. We have kept this as part of future work.

Initial guess for the optimization problem
We have described a nonlinear constrained optimization problem to find an optimal trajectory that results in a small discomfort. Because of the non-linearity and presence of both inequality and inequality constraints, it is crucial that a suitable initial guess of the trajectory be computed and provided to an optimization algorithm.
Many packages can generate their own "starting points", but a good initial guess that is within the feasible region can easily reduce the computational effort (measured by number of function and derivative evaluation steps) many times. Not only that, reliably solving a nonlinear constrained optimization problem without a good initial guess can be extremely difficult. Because of these reasons, we invest considerable mathematical and computational effort to generate a good initial guess of the trajectory.

Overview
As described in Section 2.3, a trajectory can be completely described by its length λ, orientation θ(u), and the speed v(u) for u ∈ [0, 1]. Our optimization problem is to find the scalar λ and the two functions θ and v that minimize the discomfort. We compute the initial guess of trajectory by computing λ and θ first and then computing v by solving a separate optimization problem. We emphasize that the initial guess computation process must deal with arbitrary inputs and reliably compute the initial guesses.
Before we discuss the initial guess of θ, we must discuss a genuine non-uniqueness issue. It is obvious that there exist infinitely many paths for a given pair of initial and final orientations. There exist at least two different kinds of non-uniqueness. The first kind of non-uniqueness exists because multiple numerical values of an angle correspond to a single "physical" orientation. The second kind of non-uniqueness exists because . . .   even for the same numerical values of initial and final angles, one can end up in one of multiple local minima after optimization. We now discuss these in detail.

Multiplicity of paths
Since the trajectory orientation θ is an angle, a single θ value is completely equivalent in physical space to θ ± 2nπ ∀n ∈ N. However, consider a trajectory that starts with a given angle θ 0 , and stops at orientation θ τ (where τ denotes final time). Such a trajectory will be different than a trajectory that starts off with the same orientation but stops at θ τ ±2nπ. This is because θ is continuous and cannot jump to a different value in between. Of course, the boundary condition will still be satisfied. Thus, even though the original trajectory optimization problem is specified using a single stopping orientation, we must consider multiple stopping orientations, differing by 2π, when computing the initial guess as well as solving the original discomfort minimization problem. We have called this a "parity" problem. Note that the same logic of parity applies to the starting orientation, but what matters is the difference and we have chosen to vary only the ending orientation by choosing different values of n. Figure 10 shows a few examples of this parity. It shows four paths corresponding to different n each sharing a common starting angle, but reach the destination at {−3π, −π, π, 3π}. Of course, we could create more paths by increasing the 2π difference but doing so makes the paths more convoluted and self-intersecting in general. This is because when n is too large in magnitude, θ(u) has to vary rapidly at least around some u values in [0, 1] to satisfy the larger difference in boundary conditions. For optimization purposes, we assume that the starting and ending orientations are given between [0, 2π) and we choose just three end-point orientations that give the least difference |θ τ − θ 0 |.
Apart from the multiplicity of θ curves due to parity, the optimization problem discussed Section 5.3 to compute θ as well as the full discomfort minimization problem can lead to multiple solutions even when parity remains unchanged. This occurs due to the non-linearity of constraints in Equation (15). Figure 11(a) and (b) shows two such paths A and B that start and end at the same numerical orientation but are qualitatively different. We do observe such multiple minima in practice. If we observe B more carefully, it is seen that it can be continuously deformed into B * shown in Figure 11(c) and (d). This figure clearly shows that the reason B is qualitatively different from A is because of two self-intersections − first in anti-clockwise direction and second in clockwise direction. Both "loops" cancel each others' changes in orientation. We suspect that this topological difference is the cause of multiple local minima.
Of course, this argument can be carried further and one can introduce an equal number of clockwise and anti-clockwise loops in arbitrary order and the final orientation will remain unchanged. Thus, we believe that there can be infinitely many local minima. Obviously, doing so would increase the discomfort in general and such a path will not be desirable. We try to avoid this problem by setting bounds on maximum and minimum θ when we compute the initial guess of θ. However, it is important to not ignore multiple minima if they are found within these bounds. If obstacles are present so that A is infeasible, B might be chosen even though it is longer and has more turns.
Thus, because of these two kinds of multiplicities, we use more than one initial guess when minimizing the discomfort and choose the one that has the minimum discomfort and satisfies the constraints. We discuss the details in the following section.

Initial guess of path
We compute initial guesses for λ and θ(u) using two different methods. The first method computes a θ(u) such that the trajectory has a piecewise constant curvature. This is a computationally inexpensive method and does not satisfy many of the constraints exactly. The output of this method can be used to solve the full discomfort minimization problem.
The second method computes a θ(u) and λ by solving an auxiliary (but simpler) nonlinear constrained optimization problem. Of course, now we need an initial guess for this new optimization problem! The output of the "constant curvature" method mentioned above is used as the initial guess. Unlike the first method, the output of this second method leads to trajectories that have continuous and differentiable curvature and also satisfy boundary conditions and maximum curvature constraint exactly.

Piecewise constant curvature path
In the full discomfort minimization problem, the orientation θ(u) has to satisfy the boundary conditions and Equation (15). In total, there are four constraints − two linear (those due to boundary conditions) and two non-linear (those of Equation (15)).
For computing initial guess of θ, we modify the inputs of the full optimization problem using a rotation such that initial and final position have the same y coordinate value. The initial and final orientations are also modified appropriately. Once we find an initial guess for the transformed input, it can be easily transformed back to the original configuration by the inverse rotation. This is done to allow efficient storage of precomputed θ guesses for various end-point conditions. Thus, the inputs to the initial guess generation problem are the initial and final positions, x 0 and x τ , and orientations, θ 0 and θ τ , in the rotated frame. The output will be a path length λ and a function θ(u).
We begin by choosing the value of path length λ as max(R, 2 * ∆L), where R is the minimum turning radius of the robot and ∆L = ||r τ − r 0 ||. Using this maximum takes care of the case when initial and final positions are very close to each other. In such a case, the path length is decided by the minimum turning radius constraint.
Ideally, an initial guess of θ(u) should obey the following constraints so that the constraints of the full optimization problem are satisfied: and transformed Equation (15) Consider a piecewise linear function that looks like the solid curves in Figure 12(b).
Essentially, each such curve is defined on [0, 1], is continuous, and is made of three line segments in [0, 1 3 ], [ 1 3 , 2 3 ], and [ 2 3 , 1]. The middle line segment has zero derivative. The values at 0 and 1 are known and the only variable is the function value on the middle segment. Equivalently, one can use the slope of first line segment as the variable. Let this slope be denoted by θ 1 . Then, we define θ(u) as If we use such a curve for θ(u), it will result in a circular arc, a tangent line segment, and another circular arc tangential to the middle segment, in that order. This, in turn, implies that the resulting path will have a piecewise constant curvature.
To determine θ(u), we need to determine the value of the unknown slope θ 1 . Since only one value cannot satisfy two constraints of Equation (22), we minimize to find θ 1 . Figure 12(a) shows the plot of J cc as a function of θ 1 . Depending on the boundary conditions, the shape of J cc changes but qualitatively it has the behavior as shown − oscillatory with a maximum not too far from zero. We find this maximum using a table lookup and the neighboring two minima to compute the initial guess. Figure 12(c) shows two paths using this method where θ 0 = π 2 and θ τ = π 3 . The path length is 1. As seen, the curve end-point is not too far from x−axis, and the curve satisfies the boundary condition on θ. Figure 12(d) shows the cost J cc for various constant curvature paths. The lighter shade corresponds to the minima.

Optimization approach for initial guess of path
In this second method to compute the initial guess of the path, we minimize where w := max(∆L, R), and θ must satisfy the boundary conditions, the two equality constraints of Equation (15), and the curvature constraint We do not impose the obstacle related constraints in this problem. This problem is related to the concept of "Minimum Variation Curves" Moreton (1993) which have been proposed for curve shape design. We add the curve length λ so that in the presence of multiplicities, discussed in Section 5.2 earlier, the curves with smaller lengths are preferred. This optimization problem is discretized using C 1 finite elements as described in Section 4.2, and the initial guess is the piecewise constant curvature function from Section 5.3.

Initial guess of speed
Computing the initial guess of v is relatively simpler. We solve a convex quadratic optimization problem with linear inequality box constraints to compute the initial guess. Because of convexity of the functional and the convex shape of the feasible region, this problem has a unique solution and an initial guess is not necessary to solve it. Any good quality optimization package can find the solution without an initial guess. Of course, because of the simplicity of box constraints, we can and do provide a feasible initial guess for this auxiliary problem.
First consider the case when both end-points have non-zero speed. We minimize The expressions for v (0) and v (1) come from the relation in Equation (8). The length λ is computed when the initial guess for θ is computed. Here we choose v min (u) = min(v 0 , v 1 )/2 and v max (u) is a constant that comes from the hardware limits. The function A min (u) is chosen to be the constant 10a min λ/ min(v 0 , v 1 ) where a min is the minimum allowed physical acceleration. A max (u) is chosen similarly using a max . This optimization problem is discretized using C 1 finite elements as described in Section 4.2 and leads to a convex programming problem that is easily solved. Of course, this method does not take care of cases in which one or both points have zero speed boundary conditions. If only one of the end-points has a zero speed boundary condition, we split the initial guess for v into a sum of two functions. The first one takes care of the singularity and the second takes care of the non-zero speed boundary condition on the other end-point. We now maintain only the v max constraint because v (u) is unbounded and v min = 0 naturally. If the right end-point has zero speed, we choose v(u) = v singular (u) + v non-singular (u) where v singular (u) = 16 9 2 1/3 v max u 2 (1 − u) 2/3 . This function has the correct singularity behavior and its maximum value is v max /2. The non-singular part is computed via optimization so that the sum is always less than v max . For the other case, when left end-point has zero speed, the singular part (using symmetry) is Figure 14 shows these different cases. All the imposed bounds are maintained and the initial guesses of v are smooth curves for all kinds of boundary conditions. For non-zero boundary speed, the values are 1 on the starting point and 2 on the ending point. Maximum speed is 3. Where imposed, A min = −50 and A max = 50. 24

Implementation details
We have implemented our code in C++. We use Ipopt, a robust large-scale nonlinear constrained optimization library (Wächter and Biegler, 2006), also written in C++, to solve the optimization problem to compute the initial guess of path (Section 5.3) and speed (Section 5.4). We explicitly compute gradient and Hessian problem in our code instead of letting Ipopt compute these using finite difference. This leads to greater robustness and faster convergence. We set the Ipopt parameter for relative tolerance as 10 −8 and the maximum number of iterations to 100.

Computational time for initial guess
To evaluate the reliability of our method, we constructed a set of 7500 problems with different boundary conditions and solved the full constrained optimization problem corresponding to each of the 4 initial guesses for each problem. We present an analysis of the reliability and run-time for the initial guesses for this problem set.
We generate the problem set as follows. All the problems were solved on a computer with an Intel Core i7 CPU running at 2.67 GHz, 4 GB RAM, and 4 MB L-2 cache size. Histograms of run-time for computing initial guess of speed and initial guess of path are shown in Figures 16 and 15 respectively. Histograms for all four initial guesses and all four discomfort minimization problems are shown. In all these histograms, we have removed 1% or less of cases that lie outside the range of the axis shown for better visualization. All histograms show both successful and unsuccessful cases.
For all 7500 problems, at least one initial guess was successfully computed. From Figure 15 we see that than 99% or more of initial guesses of path are computed in less than 0.2 s. From Figure 16 we see that 99% or more of initial guesses of speed are computed in less than 0.12 s.
These initial guesses were used to solve the full nonlinear optimization problem. Results for the full problem are presented in Part I. Here is a summary. Out of a set of 7500 examples with varying boundary conditions, and all dynamic constraints imposed, 3.6 solution paths, on average, were found per example. At least one solution was found for all of the problems and four solutions were found for roughly 60% of the problems. The time taken to compute the solution to the discomfort minimization problem was less than 10 seconds for all the cases, 99% of all problems were solved in less than 4 seconds, and roughly 90% were solved in less than 100 iterations.    28