Model-based control for robots has increasingly depended on optimization-based methods like Differential Dynamic Programming (DDP) and iterative LQR (iLQR). These methods can form the basis of Model-Predictive Control (MPC), which is commonly used for controlling legged robots. Computing the partial derivatives of the robot dynamics is often the most expensive part of these algorithms, regardless of whether analytical methods, Finite Difference, Automatic Differentiation (AD), or Chain-Rule accumulation is used. Since the second-order derivatives of the robot dynamics result in tensor computations, they are often ignored, leading to the use of iLQR, instead of the full second-order DDP method. In this paper, we present analytical methods to compute the second-order derivatives of Inverse and Forward Dynamics for open-chain rigid-body systems with multi-DoF joints and fixed/floating bases. An extensive comparison of accuracy and run-time performance with AD and other methods is provided, including the consideration of code-generation techniques in C/C++ to speed up the computations. For the 36 DoF ATLAS humanoid, the second-order Inverse and Forward Dynamics derivatives take$\approx 200 \mu s$, and$\approx 2.1 ms$respectively, on a 12th Gen Intel i5-12400 processor with 2.5 GHz clock-speed, resulting in a$\approx 3.2 \times$and$\approx 3.8 \times$speedup respectively over the AD approach.

Optimization-based robot control strategies often rely on first-order dynamics approximation methods, as in iLQR. Using second-order approximations of the dynamics is expensive due to the costly second-order partial derivatives of the dynamics with respect to the state and control. Current approaches for calculating these derivatives typically use automatic differentiation (AD) and chain-rule accumulation or finite-difference. In this paper, for the first time, we present analytical expressions for the second-order partial derivatives of inverse dynamics for open-chain rigid-body systems with floating base and multi-DoF joints. A new extension of spatial vector algebra is proposed that enables the analysis. A recursive algorithm with complexity of$\mathcal{O}(Nd^{2})$is also provided where N is the number of bodies and d is the depth of the kinematic tree. A comparison with AD in CasADi shows speedups of 1.5-3 x for serial kinematic trees with N > 5, and a C++ implementation shows runtimes of$\approx \mathbf{5 1} \mu \mathrm{s}$for a quadruped.

The computations of the high-order partial derivatives in a given problem are often cumbersome or not accurate. To combat such shortcomings, a new method for calculating exact high-order sensitivities using multicomplex numbers is presented. Inspired by the recent complex step method that is only valid for firstorder sensitivities, the new multicomplex approach is valid to arbitrary order. The mathematical theory behind this approach is revealed, and an efficient procedure for the automatic implementation of the method is described. Several applications are presented to validate and demonstrate the accuracy and efficiency of the algorithm. The results are compared to conventional approaches such as finite differencing, the complex step method, and two separate automatic differentiation tools. The multicomplex method performs favorably in the preliminary comparisons and is therefore expected to be useful for a variety of algorithms that exploit higher order derivatives.