Collocation methods for second and higher order systems

It is often unnoticed that the predominant way to use collocation methods is fundamentally ﬂawed when applied to optimal control in robotics. Such methods assume that the system dynamics is given by a ﬁrst order ODE, whereas robots are often governed by a second or higher order ODE involving conﬁguration variables and their time derivatives. To apply a collocation method, therefore, the usual practice is to resort to the well known procedure of casting an M th order ODE into M ﬁrst order ones. This manipulation, which in the continuous domain is perfectly valid, leads to inconsistencies when the problem is discretized. Since the conﬁguration variables and their time derivatives are approximated with polynomials of the same degree, their diﬀeren-tial dependencies cannot be fulﬁlled, and the actual dynamics is not satisﬁed, not even at the collocation points. This paper draws attention to this problem, and develops improved versions of the trapezoidal and Hermite-Simpson collocation methods that do not present these inconsistencies. In many cases, the new methods reduce the dynamic transcription error in one order of magnitude, or even more, without noticeably increasing the cost of computing the solutions.


Introduction
Direct collocation methods have proven to be powerful tools for solving optimal control problems in robotics (Posa et al., 2016;Pardo et al., 2016;Kelly, 2017;Hereid et al., 2018;Tedrake, 2023).Initially developed for aeronautics and astrodynamics applications (Hargraves and Paris, 1987;Conway and Paris, 2010), these methods have become very popular and of widespread use in the context of trajectory optimization and model predictive control, thanks to a few key advantages over indirect approaches based on the Pontryagin conditions of optimality: in general, they are easier to implement and show larger regions of convergence, and do not require estimations of the costate variables, which may be difficult to obtain accurately.Helpful tutorials and monographs like Kelly (2017) or Betts (2010), as well as open-source software for nonlinear optimization (Wächter and Biegler, 2006), numerical optimal control (Kelly, 2017;Becerra, 2010;Andersson et al., 2019), or model-based design and verification (The Drake Team, 2023), are also contributing to their rapid dissemination among the community.
Direct collocation methods involve the transcription of the continuous-time optimal control problem into a finite-dimensional nonlinear programming (NLP) problem (Kelly, 2017).The transcription is based on partitioning the time history of the control and state variables into a number of intervals delimited by knot points.The system dynamics is then discretized in each interval by imposing the differential constraints at a set of collocation points, which may coincide, or not, with the chosen knot points.The cost function is also approximated using the values taken by the variables at such points, and the NLP problem is formulated using them.Once this problem is solved, a continuous solution is built using interpolating polynomials that satisfy the dynamics equations at the collocation points.
The general formulation of most collocation methods assumes that the system dynamics is governed by a first order ODE of the form where x(t) and u(t) are the state trajectory and the control function, respectively (Tedrake, 2023).
In robotics, however, as in mechanics in general, the evolution of the system is often determined by a second order ODE of the form q(t) = g (q(t), q(t), u(t), t) , where q(t) is the configuration trajectory and q(t) is its time derivative.To apply a general collocation method, therefore, the usual procedure is to define v(t) = q(t) and write (2) as q(t) = v(t), v(t) = g(q(t), v(t), u(t), t). (3a) which, if we define x(t) = (q(t), v(t)), corresponds formally to (1).Yet, this raises a consistency issue.
Since the collocation method locally approximates q(t) and v(t) by polynomials of the same degree, imposing only at the collocation points does not grant the satisfaction of (4) over the continuous time domain.Even more striking, perhaps, is the fact that, as we demonstrate in this paper, imposing (3) at the collocation points does not imply the satisfaction of (2), not even at these points, which contributes to increase the dynamic transcription error along the obtained trajectories.This hinders the possibility to reach a correct solution since, even if u(t) produces the expected trajectory for v(t), its integration will rarely coincide with the function obtained for q(t).In other words, the state trajectory x(t) will be inconsistent in general.
In this work, we present modified versions of the trapezoidal and Hermite-Simpson collocation methods specifically addressed to solve these issues for second order systems with the dynamics in (2).The new formulations grant that the collocation polynomials fulfill the condition in (4) while satisfying (2) at the collocation points, thereby increasing the accuracy of the obtained solutions.We early presented these methods in Moreno-Martín et al. (2022), and here we also extend them to treat M th order ODEs of the form q (M ) (t) = g(q(t), q(t), ..., q (M −1) (t), u(t), t), (5) where q (M ) (t) = d M q(t)/dt M .While less common, ODEs of this kind arise more or less explicitly in flexible, elastic, or soft robots for example (De Luca and Book, 2016;Della Santina, 2020).
By means of well-established benchmark problems from the literature, we further demonstrate that the new methods reduce substantially the dynamic transcription error (in one order of magnitude or even more depending on the number of knot points) without noticeably increasing the computational time needed to solve the transcribed NLP problems.As a result, the state and control trajectories x(t) and u(t) will be mutually more consistent, thus facilitating their tracking with a feedback controller.
The rest of the paper is structured as follows.Section 2 formulates the optimal control problem to be solved and delimits the specific transcription problem that we face in this paper.To prepare the ground for later developments, Section 3 reviews the conventional trapezoidal and Hermite-Simpson collocation methods for first order systems and pinpoints their limitations on transcribing second order differential equations.Improved versions of these methods are developed in Section 4 for second order systems, and for M th order ones in Section 5.The methods are summarized and compared in Section 6, and their performance is analysed in Section 7 with the help of examples.Finally, Section 8 concludes the paper and enumerates a few points for future attention.

Problem formulation
The optimal control problem that concerns us in this paper consists of finding state and action trajectories x(t) and u(t), and a final time t f , that minimize subject to where K(x f , t f ) and L(x(t), u(t)) are terminal and running cost functions, respectively, Eq. ( 6b) is an ODE modeling the system dynamics, Eqs. ( 6c) and (6d) encompass the path and boundary constraints, x 0 = x(0), and x f = x(t f ).We note that, while Eq. ( 6b) has the appearance of a first order ODE, in robotics it often takes the form where x = (x 1 , . . ., x M ) = (q, q, . . ., q M −1 ), ( 8) so in such cases it actually encodes an M th order ODE like (5), or (2) if M = 2. Solving Problem (6) via collocation involves partitioning the time history of the control and state variables into N intervals delimited by N + 1 knot points t k , k = 0, . . ., N , then transcribing Eqs.(6a)-(6c) into appropriate discretizations expressed in terms of the values x k = x(t k ) and u k = u(t k ), and finally solving the constrained optimization problem that results.
The transcriptions of (6a) and (6c) are relatively straightforward and less relevant in the context of this paper.They can be done, for example, by approximating the integral in (6a) using some quadrature rule, and enforcing (6c) for all knot points t k .The transcription of (6b), in contrast, is substantially more involved, and will be the main subject of the rest of this paper.In particular, we seek to construct appropriate polynomial approximations of the solutions x(t) of (6b) for each interval [t k , t k+1 ].These approximations will be defined as solutions of systems of equations which, when considered together for all intervals, will form a proper transcription of (6b) over the whole time domain [0, t f ].

Methods for first order systems
Two of the most widely used transcriptions of (6b) are those of the trapezoidal and Hermite-Simpson methods, which assume no particular form for (6b).To see where these transcriptions incur in dynamical error, and ease the development of the new methods, we briefly explain how they approximate (6b) and obtain their approximation polynomials for the state.Our results match those by Betts (2010) and Kelly (2017), but we follow a derivation process that is closer to Hargraves and Paris (1987), which facilitates the transition to our new methods in Sections 4 and 5.

Trapezoidal collocation
In trapezoidal collocation, the state trajectories are approximated by quadratic polynomials.If, for each interval [t k , t k+1 ], we define τ = t−t k , we can write the polynomial approximation for a component x of the state in this interval, and its temporal derivative, as where a, b, and c are real coefficients.To facilitate the application of collocation constraints, however, we will rewrite x(t) using the three parameters ẋk+1 = ẋ(t k+1 ).( 12) Evaluating the right-hand sides of ( 10)-( 12) using ( 9) we obtain where h = t k+1 − t k , so solving for a, b, c and substituting the resulting expressions in (9a), we have Equation ( 14) is known as the interpolation polynomial, as it allows us to estimate the intermediate states for t ∈ [t k , t k+1 ], once the NLP problem has been solved.Now, following Hairer et al. (2002, page 30), we determine the three parameters of ( 14) by enforcing the initial value constraint x(t k ) = x k and two collocation constraints of the form From ( 14) we see that x(t k ) = x k by construction.As for the collocation constraints, the trapezoidal method imposes them at the knot points t k and t k+1 , so it must be where f k is a shorthand for f (x k , u k , t k ).The value x k+1 , then, is obtained by evaluating ( 14) for τ = h.This results in the constraint which ensures the continuity of the trajectory across intervals k and k + 1.Note that equations ( 15)-( 17) already form a transcription of our ODE in the interval [t k , t k+1 ] since, if x k , u k , and u k+1 were known, these equations would suffice to determine the three unknowns ẋk , ẋk+1 , and x k+1 .However, we can also substitute ( 15) and ( 16) into ( 17) to obtain the more compact expression which we recognize as the common transcription rule in trapezoidal collocation (Kelly, 2017;Betts, 2010).Observe that the continuity between the polynomials of intervals k and k + 1 is granted for the first derivative as, by construction, they both satisfy ẋk+1 = f k+1 .However, second and higher order continuity is not preserved in general.

Hermite-Simpson collocation
In Hermite-Simpson collocation, the state trajectories in each interval are approximated by cubic polynomials: By analogy with the trapezoidal method, we first express the polynomial coefficients in terms of the parameters where t c = t k + h/2, and the extra parameter ẋc is added because four parameters are needed to determine a third degree polynomial.Evaluating these identities using ( 19), solving for a, . . ., d, and substituting the expressions in (19a), we obtain the interpolation polynomial In order to determine the four parameters of (20), four conditions have to be imposed, and the Hermite-Simpson method makes this by fixing x(t k ) = x k (which holds by construction) and imposing the dynamics at the two bounding knot points and the midpoint between them: In the latter equation, f c = f (x c , u c , t c ), where x c = x(t c ), and u c = u(t c ).Moreover, the values x c that are needed in f c can be expressed in terms of the above four parameters by evaluating (20) for τ = h/2, which yields Finally, the continuity constraint between intervals k and k + 1 is obtained by evaluating (20) for τ = h: Equations ( 21)-( 25) already form a transcription of our ODE in [t k , t k+1 ], but a transcription involving less variables can be obtained by substituting ( 21)-( 23) in ( 25) and ( 24), which gives If preferred, we can also remove the dependence on f c in (26b).This is achieved by isolating f c from (26a) and substituting the result in (26b), which yields the alternative transcription Both transcriptions in ( 26) and ( 27) are called separated forms of Hermite-Simpson collocation, in the sense they both keep x c as a decision variable of the problem.They are equivalent, but the one in (27) allows us to eliminate x c by substituting (27b) in (27a), which results in a single equation that is known as the compressed form of Hermite-Simpson collocation (Kelly, 2017;Betts, 2010).While the use of a separated form tends to be better when working with a small number of intervals, the compressed form is preferable when such a number is large (Kelly, 2017).
Note that, despite the polynomial approximation into each interval between consecutive knot points is of third degree, continuity through knot points is only granted for the state trajectory and its first derivative.

Trajectory interpolation
After solving the NLP problem, values of the state and control variables at all collocation points are available.A continuous approximation to the optimal trajectory for the state is then obtained by substituting ( 15)-( 16), or ( 21)-( 23), in the corresponding interpolating polynomials ( 14) and ( 20), for the trapezoidal and Hermite-Simpson methods, respectively.The approximation of the control trajectory within each interval is obtained, in the trapezoidal case, by linear interpolation of the control values.In the Hermite-Simpson case, different options are possible.Some authors handle the midpoint control as an independent variable and use a quadratic interpolation of the three control values available in each interval (Kelly, 2017), while others prefer a linear interpolation and enforce the midpoint value to be the mean of the two bounding values (Topputo and Zhang, 2014).In this paper we follow the former option.

Downsides of the methods
In a first order dynamical system, imposing (15)-( 16) or ( 21)-( 23) grants that the system dynamics is effectively satisfied at the collocation points.The same is not true when a second order system is cast into a first order one via (3).To see why, note that the constraint q(t) = v(t) is only imposed at the collocation points, but not in between them, so that, even if the curves q(t) and v(t) coincide at such points, their derivatives may be different in them (Fig. 1).Therefore, q(t) = v(t) in general and, in particular, also at the collocation points.As a consequence, even if qk = v k and vk = g(q k , v k , u k , t k ), this does not imply that the expected relation qk = g(q k , v k , u k , t k ) is satisfied, what means that, with a transcription based on (3), the system dynamics in (2) is not granted, not even at the collocation points.This problem is solved in the second order collocation methods introduced in the next section.
< l a t e x i t s h a 1 _ b a s e 6 4 = " P r 5 a P g X p c w 4 N a U 2 P S o 9 h h y m m H U g = " > A A A B 6 H i c b V B N S 8 N A E N 3 U r 1 q / q h 6 9 L B b B U 0 l E 1 G P R i 8 c W 7 A e 0 o W y 2 k 3 b t Z h N 2 J 0 I J / Q V e P C j i 1 Z / k z X / j t s 1 B W x 8 R e 1 I n T c I J k G f y S t 6 c R + f F e X c + F q 0 F J 5 8 5 J n / g f P 4 A 4 c + M / Q = = < / l a t e x i t > t < l a t e x i t s h a 1 _ b a s e 6 4 = " t C 7 / H m e h y m k S x 2 y 0 Y M 0 9 h + h P W Z 8 / s W 2 X n w = = < / l a t e x i t > qk 6 = vk < l a t e x i t s h a 1 _ b a s e 6 4 = " W G o O 5 I Y k X h 1 9 A S V s l G g 1 t 7 H J p f c = " > A A A B 6 n i c b V B N S 8 N A E J 3 U r 1 q / q h 6 9 L B b B U 0 l E q s e i F 4 8 V 7 Q e 0 o W y 2 m 3 b p Z h N 2 J 0 I J / Q l e P C j i 1 V / k z X / j t s 1 B W x 8 M P N 6 b Y W Z e k E h h 0 H W / n c L a + s b m V n G 7 t L O 7 t 3 9 Q P j x q m T j s h a 1 _ b a s e 6 4 = " w G 9 c 0 3 7 q R s e e h y w i 5 / M l X 7 5 i Fig. 1: Inconsistencies that arise when a collocation method for first order systems is applied to a second order ODE q = g(q, q, u, t).The figure illustrates the case of the trapezoidal method, whose quadratic approximations q(t) and v(t) are depicted in blue.The red and green curves correspond to first and second derivatives of these trajectories, respectively.
A related problem of first order methods is that, when the trajectories are approximated with their interpolation polynomials q(t) and v(t), the difference v(t) − q(t) = 0 makes the state trajectory inconsistent, so that, if we try to follow it with a controller, since the configuration and velocity trajectories are incompatible, both cannot be followed at the same time.An attempt to solve this may consist in ignoring the configuration trajectory and replacing it by the integral of the velocity, but the resulting configuration trajectory may violate the problem constraints, e.g., the final configuration may be different from the expected one.Alternatively, one can try to replace the velocity trajectory by the derivative of q(t), but in this case, since the dynamic constraint satisfied at collocation point k is vk = g(q k , v k , u k , t k ), and v k = qk but vk = qk , the dynamic constraint qk = g(q k , qk , u k , t k ) will not be satisfied with the computed u k .

Methods for second order systems
To solve the inconsistency problems just explained, we propose alternative formulations for the trapezoidal and Hermite-Simpson collocation methods in which the dynamic constraints are directly imposed on the second derivative of the configuration variables, instead of on the first derivative of the state variables.By doing so, the velocity variables are not treated as independent from the configuration ones, but explicitly defined as v(t) ≡ q(t).In this way, the discrepancy between q(t) and v(t) is fully removed, and the second order dynamics is satisfied at each collocation point.

Trapezoidal method for second order systems
The essential feature characterizing trapezoidal collocation is that the dynamics is imposed just at the knot points or, otherwise said, that each interval bound is a collocation point.When the dynamics is governed by the second order ODE in (2), using the same strategy as the trapezoidal method consists in imposing (2) at each interval bound.This means that, for each interval, two constraints have to be imposed on the second derivative of the polynomial approximating each component q of the configuration.But, since the second derivative of a quadratic polynomial is constant, only one constraint could be imposed on it.This implies that the interpolating polynomial q(t) must be of degree three at least.So, we will have, for a given interval [t k , t k+1 ], To determine the coefficients a, b, c, d, we need to impose four conditions.While in the trapezoidal method three conditions were used (the value x k at the initial bound and the derivatives ẋk and ẋk+1 at the two bounds), here we will impose, in addition to the initial value q k and the second derivative at the interval bounds qk and qk+1 , the value qk of the first derivative at the initial bound.Note that, for a cubic polynomial, no more than two independent conditions can be fulfilled by its second derivative, so imposing the dynamics at the midpoint of the interval as in the Hermite-Simpson method is not possible here.Thus we will use as parameters: Evaluating these identities using (28) and solving for a, b, c, d, we can write the interpolation polynomial q(t) as: The evaluation of this polynomial and its derivative q(t) for τ = h yields and imposing the collocation constraints where g k = g(q k , qk , u k , t k ), we finally obtain the trapezoidal method for second order systems: Note that, in this case, the trapezoidal rule only applies for the velocity, but not for the configuration itself, which is given by equation (32a).
As opposed to the trapezoidal method for first order systems, the continuity between neighboring polynomials at the knot points is of second order in this case, since the collocation constraints impose the coincidence of the second derivative of q(t).Second order continuity for the configuration trajectory implies smooth velocity profiles and continuous accelerations, which are desirable properties in many robotics applications (Constantinescu and Croft, 2000;Macfarlane and Croft, 2003;Berscheid and Kröger, 2021).

Hermite-Simpson method for second order systems
Our purpose now is to impose the second order dynamics on the two bounds and the midpoint of each interval, in similarity with the conventional Hermite-Simpson method.Clearly, if we want to impose three conditions to the second derivative of a polynomial q(t), such a derivative must be quadratic at least, what implies that the polynomial must have degree four at least.Thus, we propose to approximate the configuration trajectory, and its derivatives, by Since five parameters are needed to determine the five coefficients of q(t), we will use, in addition to the three accelerations qk , qc , qk+1 , the values of the configuration coordinate q k and its derivative qk at the initial point: Solving for the coefficients a, . . ., e, we obtain the following expression for the interpolating polynomial: Evaluating ( 36) and its derivative for the value τ = h results in and imposing the collocation constraints yields where we recognize that (39b) is the Simpson quadrature for the velocity.The terms g c in these equations involve the midpoint coordinate q c = q(t c ), and the velocity qc = q(t c ), but these can be obtained by evaluating (36) and its derivative for τ = h/2, and imposing (38), which yields Note however that, since q c and qc are to be used in the evaluation of g c , we may prefer not to express them in terms of g c itself.For this we simply isolate g c from (39b) and substitute the result in ( 40) to obtain: Equations ( 39) and ( 41) together constitute a separated form of the Hermite-Simpson method for 2nd order systems.Written in this way, (41) can be replaced in the expression of g c in (39) to transcribe the problem in compressed form, which eliminates the need to treat q c and qc as decision variables.
In this collocation scheme, the continuity across knot points is also of second order due to the coincidence of the second derivative imposed by the collocation constraints, what gives rise to smooth, continuous acceleration trajectories just like in the second order trapezoidal method.

Extensions for higher order systems
Second order systems are, by far, the most common in robotics, but sometimes it may be necessary to deal with dynamical systems of a higher order M , whose dynamics is described by an ODE like (5), which we recall for convenience: q (M ) (t) = g q(t), q(t), ..., q (M −1) (t), u(t), t .(42) We next see how the new methods can be extended to transcribe (42).

The generalized trapezoidal method
To derive the trapezoidal method for M th order systems we proceed as in Section 4.1.For each time interval [t k , t k+1 ] we approximate each component q of the solution of (42) by a polynomial q(t) whose M th time derivative q (M ) (t) is linear, so its two coefficients may be determined by imposing (42) at each interval bound.This implies that q(t) must be of order M + 1.Using a 0 , . . ., a M +1 as coefficients, this polynomial, and its derivatives, may be written as . . .
or, more compactly as for j = 0, . . ., M .We then can determine a M and a M +1 by imposing the two collocation constraints where g k = g(q k , qk , ..., q The remaining coefficients a 0 , . . ., a M −1 are determined by imposing the initial value constraints for j = 0, . . ., M − 1.Using (43) we see that the left hand side of ( 48) is a j , so we readily obtain for j = 0, . . ., M − 1.Finally, by evaluating (44) for τ = h we find that the generalized versions of the Tz -2 formulas in (32) are given by q (j) for j = 0, . . ., M − 1.
One can check that, by particularizing (50) for M = 1 and M = 2, we obtain the equations of the trapezoidal method for first and second order systems given in ( 17) and (32), respectively.

The generalized Hermite-Simpson method
An analogous route can be followed to obtain a Hermite-Simpson method for M th order systems.In this case, q (M ) (t) must be quadratic in order to determine its coefficients by imposing the collocation constraints at t k , t k+1 , and t c = t k +h/2.This means that q(t) must be of degree M + 2 now, so q(t), and its derivatives, will take the form for j = 0, . . ., M .The last equation in ( 51) is and its coefficients a M , a M +1 , and a M +2 can be determined by imposing where g c = g(q c , qc , ..., q (M −1) c , u c , t c ), (56a) After simple calculations we find that As in the trapezoidal method, the remaining coefficients are determined by the initial value constraints, and we have for j = 0, . . ., M − 1.The generalized versions of Eqs.(39) can then be obtained by evaluating the expressions up to order M − 1 in (51) for τ = h,

Method Collocation equations
Tz -1 HS -1 Table 1: Collocation equations for all methods of the trapezoidal and Hermite-Simpson families.For the Tz -M and HS -M methods, we provide the general equation of q (M −l)

k+1
(where l is meant to run up to M ) but also the particular instances of this equation for l = 1, 2. The equation of q (M −l) c , and its instances for l = 1, 2, are also provided in the HS -M method (where, again, l = 1, . . ., M ).This arrangement allows us to realize that, within each family, the equations for the same l coincide for all orders.

Family
Table 2: Number of variables (nv), equations (ne), and degrees of freedom (n DOF ) in the two families of methods.
As it happens in the Hermite-Simpson method for 2nd order systems, g c in (57) requires the midpoint values q c , qc , . . ., q , but these are easily obtained by evaluating (51) for τ = h/2, which results in The terms a M +1 and a M +2 in (60) involve g c and thus the midpoint coordinate q c and its derivatives.However, we can remove the dependence of q (j) c on g c by using the last equation in (59), which is By isolating g c from this equation we have which we can substitute in the expressions of a M +1 and a M +2 involved in (60).With these substitutions applied, ( 59) and ( 60) form a separated form of the Hermite-Simpson method for M th order systems.The condensed form is finally achieved by substituting the new version of (60) in the expressions of a M +1 and a M +2 in ( 59).Again, one can verify that, for M = 1 and M = 2, Eqs. ( 59) and ( 60) yield the Hermite-Simpson formulas for first and second order systems given in (26), and in ( 39) and ( 41), respectively.

Comparison of the methods
Table 1 summarizes the equations for all methods of the trapezoidal and Hermite-Simpson families.For short, we refer to the methods in each family by Tz -and HS -, followed by a number that indicates the order assumed for the system dynamics.For the general Tz -M and HS -M methods, the table provides the equations for q (M −l) k+1 , as well as q (M −l) k+1 and q (M −l) c , respectively, where l runs from 1 to M in all cases.We also specialize these equations for l = 1, 2, so the reader can realize that, within each family, the equations for a same value of l coincide for all orders.
In the table, the equations for the Hermite-Simpson methods are given in their separated form, and in the HS -M method we show those that result from applying the manipulations described in Section 5.2.

Problem size
It is not difficult to see that, for all methods in a same family, the number of variables (n v ), equations (n e ), and degrees of freedom (n DOF ) is the same in the resulting transcriptions of Problem (6).If n x and n u are the dimensions of x and u, and n b is the number of boundary constraints in Eq. (6d), we obtain the values in Table 2.Note that for an M -th order ODE, the state x includes the configuration vector q and its derivatives, so that n x = M n q , where n q is the dimension of q.The improved formulas, as compared to those of the Tz -1 and HS -1 methods, neither increase the problem size, nor reduce the freedom to find the optimal solution.Moreover, since the dynamic function must be evaluated at each collocation point, the number of evaluations is the same in all methods of a same family, so the new methods should not increase the cost of each iteration when solving the transcribed NLP problem.This point is also supported by the computational experiments that we present in Section 7.

Accuracy of the approximations
While the new methods introduced in this paper are explicitly designed to preserve the consistency between the configuration trajectory and its derivatives, a further question is how the application of these methods may affect the accuracy of the solution approximations and its rate of convergence as h → 0. To answer this question, we draw upon the concept of order of accuracy (Betts, 2010), or simply order (Hairer et al., 2002) of a collocation method, which, in turn, relies on the definition of local error of an approximation.
The local error k of a collocation method at interval k is defined as the difference between the computed value q k+1 and the value for t = t k+1 of the exact solution of the ODE, q(t), that passes through the computed point q k .If a collocation method approximates the solution with polynomials of degree d, we have for interval k: and the computed value q k+1 is obtained by setting t = t k + h, which means setting τ = h: On the other hand, the Taylor expansion of the solution q(t) that passes through the computed point q k is q and evaluating for τ = h we have: thus, the local error k is given by the difference of the two Taylor expansions ( 63) and ( 64): A collocation method is said to have order of accuracy p if the sum of the first p terms of ( 65) is zero.Note that this does not imply that each term vanishes by itself: when h takes a specific numerical value, different non-null terms of the sum may add to zero.In the hypothetical case that the exact solution q(t k + t) was a polynomial of degree p, a method of order p would have no local error.For this reason, the order of accuracy is also called the degree of exactness (Dahlquist and Björck, 2008), and an equivalent definition for it is that a collocation method has order of accuracy p if it is exact for all polynomials of degree ≤ p.
In the limit, when h → 0, the error of the approximation will converge to zero.The rate of this convergence is an important property of a method, and is directly given by its order.If a method has order p, the lower power of h appearing in k is p + 1, so that, when h → 0, the local error decreases as h p+1 : In all collocation methods discussed here, the interpolating polynomial used in each interval has degree d = M + s − 1, where M is the order of the ODE and s is the number of collocation points of each interval.This value d ensures the unique determination of the d + 1 polynomial coefficients given the s collocation constrains and the M initial conditions.In the event that the exact solution happens to be a polynomial q(t) of degree d, it must necessarily coincide with the interpolatory polynomial, and q k+1 will coincide with the exact value q(t k + h).This shows that the order of accuracy of any method is at least p = d = M + s − 1, so the orders of Tz -1 and HS -1 are at least 2 and 3, respectively, while the orders of Tz -2 and HS -2 are at least 3 and 4.However, these lower bounds can be surpassed in some cases.For example, the HS -1 method is known to have order 4 (Hairer et al., 2002), while its corresponding lower bound is 3.This is because it takes advantage of a special property of a family of polynomials of fourth degree.It can be proved that any fourth degree polynomial satisfying takes always the same value q(t k + h) = q k+1 .Since the only third degree polynomial satisfying these same conditions is a particular case of this family, it satisfies q(t k + h) = q k+1 , so its order of accuracy is 4.Even if the HS -2 method does not benefit from a similar property, it is granted that its order is at least as large as that of HS -1, i.e., 4.
So, we can say that the order of accuracy of the presented methods for second order systems is equal or higher than that of the corresponding methods for first order systems.In general, for the same number of collocation points, a method for M th order systems has this lower bound M − 1 units higher than the corresponding method for first order systems.

Test cases
The performance of all methods is next evaluated and compared using three trajectory optimization problems shown in Fig. 2. We refer to them as the cart-pole, bipedal walking, and ball throwing problems, respectively.The first two problems are solved and documented in detail by Kelly (2017), and thus serve to compare our results with those published in the literature.The third problem is proposed by the authors to illustrate the methods on a widely-used robot with a complex dynamics.Since analytical solutions for these problems are not available, we compare the different methods by computing the dynamic transcription errors they produce.To this end, we define the following errors relative to the dynamics constraints.
The first order dynamic error for the q i coordinate is defined as In general, this error is non-null in Tz -1 and HS -1, as these methods do not enforce v i (t) = qi (t) for all t.For the same coordinate, the second order dynamic error is We found this error to be more meaningful than the ε [1] qi (t) error reported by Kelly (2017), since it reflects the deviation from the actual system dynamics, which is expected to be minimized with the optimization process.Another error that is useful to define when all coordinates of q have the same units is the joint error for r = 1, 2. Finally, to summarize the error functions in just one number, we can compute their integrals over [0, t f ]: To compare the methods, we have implemented them in Python, using the toolbox CasADi (Andersson et al., 2019) to solve the constrained optimization problems that result.CasADi provides the necessary means to formulate such problems and to compute the gradients and Hessians of the transcribed equations using automatic differentiation.These are necessary to solve the optimization problems, a task for which we rely on the interior-point solver IPOPT (Wächter and Biegler, 2006) in conjunction with the linear solver MUMPS (Amestoy et al., 2001).The whole implementation can be downloaded from https://github.com/AunSiro/optibot, but the reader can also reproduce the results for the cart-pole and bipedal walking problems through interactive Jupyter notebooks online (Moreno-Martín, 2022a,b).The execution times we report have been obtained on a single-thread implementation running on an iMac computer with an Intel i7, 8-core 10th generation processor at 3.8 GHz.

The cart-pole swing-up problem
The cart-pole system comprises a cart that travels along a horizontal track and a pendulum that hangs freely from the cart.A motor drives the cart forward and backward along the track.Starting with the pendulum hanging below the cart at rest at a given position, the goal is to reach a final configuration in a given time t f , with the pendulum stabilized at a point of inverted balance and the cart staying at rest at a distance d from the initial Fig. 2: Benchmark problems.Left: A cart-pole system that has to perform a swing-up motion.Center: a walking biped whose periodic gait must be optimized (the three snapshots illustrate the motion that occurs between the toe off and heel strike events defining a period of the gait).Right: A 7R Panda robot that has to pick a ball at the shown configuration, and throw it from the same configuration at 10m/s horizontally.
position.The cost to be minimized is where u is the force applied to the cart, and we adopt the same dynamic equations and problem parameters as in Kelly (2017).An animation of the solution obtained with HS -2 and N = 25 can be seen in https://youtu.be/M0ivg8s-I8.Figure 3 compares the errors ε [1] qi (t) and ε [2] qi (t) obtained by the methods for the variables q 1 and q 2 shown in Fig. 2. The number N of intervals used in the comparison is 50 for the trapezoidal scheme, and 25 for the Hermite-Simpson one.This yields a fair comparison, as then the number of collocation points, variables, and degrees of freedom of the optimization are the same in all NLP problems (cf.Table 2).The plots corresponding to the first order error ε [1] qi (t) of Fig. 3, in the first and third row respectively, confirm that Tz -1 and HS -1 present a non-negligible first order error, while in Tz -2 and HS -2 this error is exactly zero as expected.
The plots in the second and fourth rows clearly show a discontinuity at the knot points of the second order error ε qi (t) for Tz -1 and HS -1, reflecting the discontinuity of q(t) at these points.In contrast, for Tz -2 and HS -2, the error functions are continuous and vanish at the collocation points, evidencing that, as anticipated in Section 3.4, the system dynamics is exactly satisfied at all collocation points for the new methods, but not for the conventional ones.
The figure also shows the dramatic reductions of ε [2] qi (t) for the new methods when compared with the corresponding Tz -1 and HS -1 ones.The numerical evaluation of the results appears in Table 3, which provides the computation times t c and the integral errors E [2] qi for this problem.It can be seen that the errors E [2] qi are almost one order of magnitude lower for Tz -2 and HS -2 than for their counterparts Tz -1 and HS -1, despite using a very similar computation time.It is interesting to see that the errors qi achieved by Tz -2 are about a half of those of HS -1 for the same number of collocation points.The comparison is relevant since both methods use polynomials of the same degree to approximate q i (t).
All joints are powered by torque motors, with the exception of the ankle joint, which is passive.Like the cart pole system, therefore, this robot is underactuated, but it is substantially more complex.
For this example we use the dynamic model given by Kelly (2017), which matches the one in Westervelt et al. (2003) with parameters corresponding to the RABBIT prototype (Chevallereau et al., 2003).We assume the robot is left-right symmetric, so we can search for a periodic gait using a single step, as opposed to a stride, which involves two steps.This means that the state and torque trajectories will be the same on each successive step.
As in Kelly (2017), we define q as the vector that contains the absolute angles of all links relative to ground, while u encompasses all motor torques.Also as in Kelly (2017), and similarly to the cart-pole problem, our goal is to find state and action trajectories x(t) and u(t) that define an optimal gait under the cost Several constraints are added to ensure a feasible gait.First of all we require the gait to be periodic, so where x 0 and x f are the initial and final states of the robot, and f H is the heel-strike map.The states x 0 and x f are unknown a priori, but constrained by (74), which is the particular form of the boundary constraint (6d) in this case.To construct f H it is assumed that, at heel strike, an impulsive collision occurs that changes the joint velocities but not their angles, and that, as soon as the leading foot impacts the ground, the trailing foot loses contact with it.The collision conserves angular momentum but introduces an instantaneous drop of kinetic energy in the system (Kelly, 2017).Next, we require the robot to march at a certain speed, which is achieved by setting the final time of the period to t f = 0.7s, and the length D in Fig. 2 to 0.5m.We also constrain the vertical velocity component of the trailing foot to be positive at t = 0, and negative when it touches the ground for t = t f .Finally, we require the swing foot to be above the ground at all times.An animation of the solution we obtain can be seen in https://youtu.be/dtS-WbESiW0.
Figure 4 shows the second order errors ε [2] for the different collocation methods.As before, the number of intervals used in the trapezoidal cases is twice those used in Hermite-Simpson ones so as to have identical number of collocation points and have balanced comparisons.The results are qualitatively similar to those of the cart-pole, though here the error diminution obtained by the new methods is even more accentuated.As we can see in Table 4, the integral second order error E [2] of HS -2 improves in more than one order of magnitude that of HS -1 even using a slightly lower computation time.In the case of Tz -2, its improvement over Tz -1 is still higher, reaching a reduction factor near 66, and using a computation time only slightly longer.

A ball throwing problem
As a final example, we apply the methods to compute an object-throwing trajectory for a 7R Panda manipulator.The robot is initially at rest, grasping a ball with its gripper, and its task is to throw the ball from the same configuration after 1 second, with an horizontal velocity of 10m/s.Since the dynamic model is complex in this case, we rely on the advanced dynamics engine Pinocchio (Carpentier et al., 2019) to compute the f k and g k values in the collocation formulas.This engine implements the forward dynamics algorithms by Featherstone (2008) in C++, which speeds up the computations considerably.As for the cost function, we use where K a is a small value that we fixed to 0.1 in our runs.While the first term in the integrand penalizes large control torques, the second helps to achieve smoother trajectories for the state.
To compare the methods on an equal footing, in all runs we feed the NLP solver with an initial guess that allows the convergence to a similar solution.This guess is obtained using the Tz -1 method with N = 25, initialized with u k = 0 for all k, and using x k values that interpolate the initial state, a guessed state for t = 0.5s, and the final state.The trajectory obtained via Tz -1 is then used to warm start all methods in the comparisons.As a reference, Fig. 5 shows the trajectory obtained using HS -2 and N = 100.Note how the robot performs a circular motion, exploiting gravity to gain momentum, so as to get back to the launch point with the required speed.
Fig. 6 compares the dynamic error ε [2] for the trapezoidal and Hermite-Simpson methods (left and right plots, respectively).As in the previous problems, the new methods notably outperform Ball throwing: Dynamic error [2] , N = 50 HS-2 HS-1 knot & collocation points collocation points Fig. 6: Second order dynamic errors for the ball throwing problem.
the conventional ones in terms of this error.The integral errors E [2] corresponding to these figures can be seen in Table 5, together with those of E [1]  and t c , confirming similar trends as in the earlier problems.

Performance scaling with N
To evaluate the performance of the methods when the number N of intervals increases, a series of experiments have been conducted by progressively rising N from 20 to 200.Each experiment has been launched several times and the average of the integral second order errors and computation times are represented in Fig. 7.For the bipedal walking and ball throwing problems we provide E [2] .For the cart-pole problem we cannot use E [2] as its coordinates q 1 and q 2 have different units.Instead we provide only the plot of E q1 , since the one of E [2] q2 is very similar.
In the three test problems, the best results for the second order error (shown on the left hand side of Fig. 7 in logarithmic scale) are those of HS -2, which, in many cases, improve the results of HS -1 in about one order of magnitude, or even more, and the improvement rate tends to increase with the number N of intervals.The same behavior is observed for Tz -2 with respect to Tz -1.Interestingly, in all cases the performance of Tz -2 produces, for the same number of intervals N , only about twice the error of HS -1, and this rate is kept rather constant with N .However, a more balanced comparison would be to look at experiments with equal number of collocation points, what means to compare each N value of HS -1 with the 2N value of Tz -2.A close look at the plots will convince the reader that this comparison gives equal or better results for Tz -2 in all cases.
The right hand side plots in Fig. 7 show the growth of the computation times with the number N of intervals.The plots consistently show that the difference in computation time between a method for first order systems and the corresponding method for second order systems is not relevant.In all cases, the growth is nearly linear in N , but the increase rate is higher for the HSmethods than for the Tz -ones.Despite the different complexity of the three problems analyzed, reflected in the different time scales involved, in all cases, the increase rate of the HS -methods is nearly twice that of the Tz -methods.In other words, the increasing rate is very similar for all methods when comparing the computation times for the same number of collocation points.

Conclusions
Trapezoidal and Hermite-Simpson collocation methods are very popular in the robotics community.However, they are conceived for dynamical systems of first order, while the dynamics of the systems found in robotics are often M th order, with M > 1.The transcription of an M th order ODE as a first order one has the unexpected effect that the dynamic equations are not actually imposed at the collocation points.Properly imposing the M th order constraints at the same such points as in the original algorithms requires increasing the degree of the polynomials approximating the configuration trajectory, while keeping the implied degrees for its time derivatives.This is achieved with the methods we propose, which grant the functional consistency between the trajectories of all the state coordinates, not only at the collocation points, but also along the computed trajectory.Using benchmark problems of increasing complexity, we have also shown that the new methods provide trajectories with a much smaller dynamic error than those of conventional methods, despite they require a comparable amount of computation time.This implies that the obtained trajectories will be more compliant with the system dynamics, so they should be easier to track with a feedback controller.Moreover, the trajectories of the new methods are M times differentiable, so in addition to enjoying smooth velocities, their accelerations will be continuous, or even the jerk if M ≥ 3, which are very desirable properties from a control perspective.
Points that deserve further attention are the extension of these ideas to pseudospectral collocation methods, which we initially explored for M = 2 in Moreno-Martín et al. (2022), or generalizations to deal with constrained multibody systems (Posa et al., 2016;Bordalba et al., 2023), floating-base systems (Tedrake, 2023, Chapter 17), or limits on jerk or higher-order derivatives of the trajectories.

Fig. 3 :
Fig. 3: Cart-pole problem: Plots of the first and second order dynamic errors ε

Fig. 5 :
Fig.5: Trajectory obtained for the ball throwing task.The robot, initially at rest, progressively gains momentum assisted by gravity, so as to get back to the initial configuration to throw the ball at the required speed.An animation of the trajectory can be seen in https://youtu.be/NsEv6JrSN8c.

Fig. 7 :
Fig. 7: Dynamic error (left) and optimization time (right) for the three test problems, as N is increased.

Table 3 :
Performance data for the cart-pole problem.

Table 4 :
Performance data for the bipedal walking problem (E

Table 5 :
Performance data for ball throwing problem (E