# Response To Reviews

Accelerating Second-Order Differential Dynamic Programming for Rigid-Body
Systems
John N. Nganga and Patrick M. Wensing

We would like to thank the editors and reviewers for taking the time to
carefully review our manuscript and for all of the helpful feedback. The
feedback has benefited the paper toward refining our contributions and
improving clarity throughout.

Detailed responses to reviewers are provided in this response document, with
changes in the manuscript highlighted in blue.

##  1 Associate Editor

he paper has been evaluated by a good number of Reviewers, who all agree on
the fact that the work is very good/excellent. The AE concurs with the
Reviewers and thinks that the work may be material for a fruitful discussion
at IROS and a good publication in the IEEE RA-L. However, the AE also believes
that the Authors should revise and resubmit the manuscript, after addressing
all Reviewers’ comments and clarifying the following. Since the computational
cost is evaluated in terms of the number of degrees of freedom, it is not
clear if this solution is actually advantageous, especially for traditional
robots which have fewer DoFs (e.g. 7). As some Reviewers observed, multiplying
constants also counts in the complexity, especially for small numbers.

We thank the Associate Editor for their comments in regards to this work, and
for recruiting such constructive reviewers. Longer term with this work, we are
highly interested in legged robots and we aim to apply our work to quadrupeds.
This paper represents a fundamental theoretical undertaking toward this aim.
In a whole-body MPC setting, legged robots are of high state dimensions
$(>=7)$. For example, [tassa2012synthesis] and [koenemann2015whole] use DDP
for a $22-$ and $25-$DoF humanoid respectively whereas [li2020hybrid] use a
simplified $7-$DoF quadruped in 2D. A complete model of a quadruped has 18
DoFs [MiniCheetahThesis].

While we agree with the reviewers that the multiplying constants also matter
in addition to the $\mathcal{O}(\cdot)$ complexity, we note that the highest
order costs, as captured by the $\mathcal{O}(\cdot)$ complexity, represent the
heaviest computational cost when scaling the algorithm to complex systems. In
this work, we mainly aim to reduce this heaviest computational cost. In fact,
the total computational complexity of derivative computations in DDP is
typically $\mathcal{O}(Nn^{3})$ where by $N$ is the horizon length and $n$ is
the DoFs. The main technical contribution of this work is the reduction of
this complexity to $\mathcal{O}(Nn^{2})$. While it is possible that an
$\mathcal{O}(Nn^{3})$ algorithm may be cheaper than a $\mathcal{O}(Nn^{2})$
algorithm for small $n$, with our updated presentation, we show (Fig. 4) that
that is not the case for us. Our code previously had extra overhead in places
from the interface between MATLAB and CasADi, and this overhead was skewing
some of the results for small $n$ previously. We have taken care in the
revision process to re-run the results to ensure consistency across them, and
fairness across the comparisons. We thank the AE and reviewers for the
feedback that led us to correcting our minor mistake, since correcting it
allows our contributions to be shown more clearly.

nother index to be considered is the basin of convergence for the algorithm: a
slower convergence but with guarantees from further initial guessed is
sometimes preferable. This should be explained and addressed in the revision
especially in view of traditional robots.

Towards this aim, we randomize the initial control sequence for the $7-$link
KUKA LBR manipulator and present the following findings along with Fig. 1

[2.5em]0em  Since the solution of an OCP is dependent on the initial
conditions, we randomize the initial controller using Ornstein-Uhlenbeck
process noise for the $7-$link KUKA LBR manipulator and solved the OCP problem
using either iLQR/DDP in $40$ separate instances. The first pane of Fig.
LABEL:fig:RandomCtrl illustrates the progression of that manipulator from an
initial configuration to the balanced upright configuration along an optimal
trajectory. Subplot (a) of Fig. LABEL:fig:RandomCtrl illustrates the empirical
probability density function (pdf) of the number of iterations over those
optimizations. Subplot (b) illustrates the pdf of the log of the final cost of
each optimal solution. As illustrated in subplot (a), in $95\%$ of the
simulations, iLQR had a higher number of iterations. On average, iLQR had
three times as many iterations as DDP. From subplot (b), we note that most of
the simulations regardless of DDP/iLQR converged to similar solution; in fact,
DDP converged to a different solution than iLQR in only three instances.
Figure LABEL:fig:RandomCtrl illustrates that the inclusion of second-order
information in DDP will result in similar converged solution and in fewer
iterations. This result is powerful in that the addition of second-order
information results in algorithm that converges in less iterations and to the
same minima.

Figure 1: Simulations on $7-$link KUKA LBR manipulator considering randomized
initial controller given by an Ornstein-Uhlenbeck process noise.

We also note that the DDP/iLQR iteration process can be viewed as a complex
discrete-time nonlinear process mapping a control guess at one iteration to a
control guess at the next. While theory points toward DDP converging more
rapidly to a local optima, we wish to note that the positive empirical results
herein are likely system dependent, i.e., this results may not hold for
another system, or the benefits could be more striking.

oreover, a suggestion that the AE wants to highlight, in favor of the
approach, is that recent works have shown the equivalence between soft-bodied
robots and rigid ones with higher numbers of joints/degrees of freedom. This
has already allowed the application of some classical control methods to soft-
bodied robots. The Authors may consider to investigate whether or not in these
settings the advantage of their approach compared to e.g. traditional iLQR is
more apparent

We thank the Associate Editor for their analysis towards this open research
question. Multi-segment soft-body robots can indeed be matched to an augmented
rigid-body model (c.f., the excellent work in [softToRigid]). If there are
other references that the AE or reviewers would suggest, we would be happy to
include them, as this is not our main area of expertise. While the direct
applicability of our methods to such augmented models is clear, this topic
would be worth its own research analysis. For example, the augmented rigid-
body model for soft robots introduces elastic and dissipative terms that
directly encode the impedance of the soft robot [softToRigid]. While these
effects could be incorporated as external forces into the model, it is unclear
how that would impact iLQR vs. DDP derivative calculation time, and that would
be interesting to benchmark more fully. It would also be quite interesting to
study how generalizations of the RNEA and ABA to soft robots
[DiscreteCosseratSoftRobots] could be used to go yet beyond the rigid-body
assumption in this current work. We have updated the conclusion to highlight
these exciting future prospects for the dynamic control of soft robots.

[2.5em]0em Finally, our work may find applicability when working to control
soft robots. Many multi-segment soft-body robots can be modeled assuming
piece-wise constant curvature (PCC) [softToRigid], approximating PCC with a
high-DoF augmented rigid model [softToRigid], or considering discrete Cosserat
models [DiscreteCosseratSoftRobots]. The ideas herein may find applicability
for these models due to their dynamics equations taking a similar structural
form (e.g., (LABEL:eq:id)) as rigid-body systems. Augmented rigid models are
most directly applicable in 2D, but generalizations of the RNEA and ABA
[DiscreteCosseratSoftRobots] for the 3D continuum case also present
interesting future avenues to broaden the application of this work.

nother important issue relates to how general is the proposed solution. Indeed
by clarifying this, authors may actually highlight better the advantages of
the method.

We note that the first contribution whereby we recognize the structure of the
dynamics equations (and partials thereof) and use reverse-mode to calculate
the needed quantities is general. But the second contribution where we extend
the relation of the first-order sensitivity of $\operatorname{ID}$ to
$\operatorname{FD}$ to the second-order case is specific to rigid-body
systems. With this clarification in mind, we edit the contribution as follows:
[2.5em]0em  This work presents and combines several advances for reducing the
computational complexity of computing second-order dynamics sensitivity
information in DDP. The final result is a method for computing this
information with the same computational complexity as first-order dynamics
derivatives (e.g., as in iLQR). The contributions are (I) the use of reverse-
mode automatic differentiation (AD) [grievank2000principles] to compute
second-order derivatives needed in DDP. This contribution is general to
discrete-time dynamic systems and enables computation reductions compared to
methods that explicitly evaluate a derivative tensor for the dynamics.
Further, we show (II) how second-order information related to the forward
dynamics can be related to associated information from the inverse dynamics,
akin to first-order results in [carpentier2018analytical,
jain1993linearization]; this contribution is specific to rigid-body models.
Lastly, we (III) introduce a modification to the Recursive-Newton-Euler
Algorithm (RNEA) that supports this process and further reduces computational
demands. Figure LABEL:fig:Fig1 overviews a benchmark of the proposed methods
against iLQR and against DDP approaches that explicitly compute derivative
tensors.

We have also edited the title to ‘Accelerating Second-Order Differential
Dynamic Programming for Rigid-Body Systems’ to clarify the overall method we
propose is specific to rigid-body systems.

Final Remarks:

We humbly acknowledge that responding to the many helpful reviewers has
required judicious consideration for how to spend the limited revision time of
RA-L while also considering its page limitations. However, we have done our
best to provide full answers in this document in cases where space did not
allow expansion in the main text.

As a final note, we would graciously accept further feedback from the
reviewers as to whether referencing Figures 2 and 3 early in the development
(to provide a roadmap for the reading) was helpful, of if referencing them
early (before 100% of their details were presented in the text) was more
confusing.

##  1 Reviewer 3

n this work the authors propose an efficient approach for computing higher-
order dynamics derivative terms for use in the differential dynamic
programming algorithm. This derivative information can be utilized/acquired
without explicitly forming the associated third-rank tensors. As a result, the
authors demonstrate that differential dynamic programming (DDP) with complete
second-order information can be performed at a computational cost comparable
to iterative LQR (iLQR), which only used first derivatives of the dynamics.
This approach is demonstrated on a number of n-link pendulum systems.
Strengths: Being able to perform DDP at the same computational cost as
iterative LQR is a novel contribution to the motion planning community. Near a
local minimum, DDP exhibits quadratic convergence compared to iLQR’s
superlinear convergence. The paper is well written and should be easy to
understand for a researcher working on motion planning or control. Weakness:
The paper can be greatly strengthened by addressing the following items.

We thank the reviewer for their comments and improve upon the document as
follows.

n Fig 1., what is ‘required time’? I assume this is the time for the algorithm
to converge. What is the termination criteria?

In Fig. 1, the ‘required time’ refers to the time to solve either iLQR/DDP to
convergence whereby the same convergence criteria was used for all algorithms.
The convergence criteria was indicated by a negligible ($<10^{-9}$) reduction
in the cost function between subsequent iterations. We kindly refer the
reviewer to the opening of Section IV for further detail.

hat is the initial guess provided to the solver? Without these specifications,
this plot doesn’t provide the reader much insight into how these methods
actually compare

In Figure 1 and Figure 5, we used a dissipative “swing-down” controller as the
initial control trajectory for the swing-up optimization starting from a fixed
initial condition near the downward configuration. We added a note about the
dissipative initial controller to the caption of Figure 1.

In Figure 6 and Figure 8, we used a dissipative controller as the initial
control trajectory and randomize the initial state vector around the ‘down’
position of the pendubot.

nother weakness is that the paper focuses exclusively on n-link manipulators
(that I assume are fully actuated). The paper should consider additional
canonical systems, e.g., cartpole, acrobot, quadrotor. From the n-link
pendubot examples alone it is not possible to determine if this tensor-free
approach works more generally.

Due to page limit constraints, we chose to focus on underactuated $n-$link
pendubots whereby the last link is unactuated. For the $7-$ link KUKA LBR
manipulator example, it is indeed fully actuatuted. These canonical systems
capture complex dynamics and nonlinearity aspects well, and the general
algorithmic contributions generalize to other rigid-body systems as well. With
so many reviews to juggle, we are sorry that we could not run dedicated
results for the review response to show this generality, even though there is
not space in the paper. However, we would like to emphasize that the
theoretical complexity results for the derivatives calculations (our main
contribution) will hold for these other systems.

nother comparison the authors should provide is how well the tensor-free
approach compares to iLQR when the current iterate/solution is far from the
minimum–this probably is related to situations where regularization is
required. During trajectory generation for this case, my intuition is that
iLQR will be able to more rapidly find a coarse solution compared to DDP. I
imagine that tensor-free DDP would be great for model-predictive control.

In Section LABEL:sect:Results_control, we note that we randomize the initial
state vector around the ‘down’ position of the pendubot. This choice allows
DDP/iLQR to explore the state space more fully. Our results show that DDP and
iLQR most often converge to the same solution but that DDP does it in fewer
iterations. As the reviewer suggests, iLQR is able to make more progress in
the initial iterations, as shown in Fig. LABEL:fig:SuboptimalityIteration of
the manuscript, since it is able to avoid regularization. However, this may
not hold for other more general non-convex cost functions.

We agree that DDP will be great for model-predictive control, and indeed this
is our next step. We look forward to sharing those results in due course.

In terms of comparing iLQR vs. DDP for guesses that are far off from optimal,
we chose to prioritize the analysis suggested by the AE (see, AE comment
AE.2). We show that DDP and iLQR usually converge to the same solution, even
when the initial guess is far away.

verall, the paper presents a novel contribution for motion planning using
differential dynamic programming. However, the paper can be greatly
strengthened by improving the results section with analysis of additional
systems and addressing scenarios where DDP may perform poorly compared to
iLQR.

We thank the reviewer for their positive comments. We have broadened the
results section to include more analysis on the $7-$link KUKA LBR manipulator
in LABEL:sect:Results_derivs including analysis that quantifies the quality of
the final solution. For this case, while DDP has fewer iterations, there are
some cases where it has many more iterations and gets stuck in a poor local
optimum. It is quite possible that this behavior is system dependent, but we
do agree that there are some cases (for some initial guesses) where iLQR may
be preferable. Regardless, our new tensor-free DDP is always preferable to the
conventional DDP approach.

##  4 Reviewer 6

he paper describes a method for optimizing second order derivative
computations required for Differential Dynamic Programming applied to rigid
body dynamics. By leveraging the structure of DDP applied to rigid body
dynamics the computational complexity is reduced compared to a naive
implementation of the second order approximation while maintaining the
quadratic convergence properties of DDP. Overall the paper is very clearly
written, and the approach and the evaluation are sound. The authors show that
the calculation of the second order approximation is indeed optimized,
reducing the complexity. The results also show that the convergence properties
of DDP are maintained. Improvements can be made to the title ”Tensor-Free
Second-Order Differential Dynamic Programming”. The title is suboptimal for
two reasons: 1\. The method is only applicable to Rigid Body Dynamics and not
general systems. This should definitely be included in the title, especially
for cross-readers from other fields. 2\. ”Tensor-free” is also not clear
without context. E.g. iLQR has several matrices indexed over time, which can
be represented as tensors, depending on the implementation. So ”tensor free”
in the title sounds like an implementation detail but in fact it is not, it is
an optimization of the computation of the second order approximation.

We thank the reviewer for their comments and we agree that the previous title
lacked this specificity. We re-title the paper as ‘Accelerating Second-Order
Differential Dynamic Programming for Rigid-Body Systems’ to better capture our
contributions.

he results show that there is an advantage to the proposed method over iLQR
but it seems that it is limited. One reason for that could be that rigid-body
dynamics systems can locally be well approximated linearly. Given that the
method is limited to these systems, a problem that is ”more nonlinear” is hard
to find. Maybe a case with contacts could help here?

While we nominally agree with the reviewer, we kindly acknowledge that it is
hard to rigorously quantity ‘nonlinearity’ and therefore we used canonical
models to form a comparative basis. Even so, for these models, Fig.
LABEL:fig:RandomCtrl shows that DDP will converge in about half the iterations
of iLQR.

The work as presented herein is for models without contacts and we aim to
extend this work for models with contacts. We are close on the math, but the
generalization to this broader class is quite non-trivial. You may find a nice
description for the contact treatment for the iLQR case in the thesis of Rohan
Budhiraja [budhiraja2019multi]. Overall, we appreciate (and share) the
optimism of the reviewer regarding the use of our methods in the contact case,
and hope to be able to report on further results soon.

he comparison to analytic derivatives (as e.g. implemented by Carpentier et
al.) would be interesting and this is mentioned as future work. Understandably
this is not included in the paper yet but would have made it the contribution
much stronger.

We have indeed extended our work to accumulate the second-order derivatives in
similar fashion as the first order ones in [carpentier2018analytical], but
this has only been in the past week :) The final form of the second-order
partials of modified RNEA is identical in structure (but quite different in
details) as the first-order partials of regular RNEA. [The modified RNEA
second-order partials are more complicated, since they are second order, but
computing them is only a constant factor more costly than the first-order RNEA
ones.] While the extension to contacts is non-trivial, the second-order
derivatives are even more so. We will be happy to share these algorithms open-
source with the community (via planned additions to the phenomenal Pinocchio
software package) after we have them fully tested.

nother small point is that the authors use Euler forward integration, which,
for derivation purposes, is fine. However, given that it is not energy
conserving and not really state of the art anymore, it would be great for the
results to use a semi-implicit integration scheme. Given that there are
toolboxes such as Acado that allow for C code generation of state-of-the-art
optimal control solvers based on systems implemented in CasADi, a wider
comparison would be great.

This work is straightforwardly applicable to other explicit integration
schemes and our choice of Euler forward integration is merely for a simple
first implementation sake and not because of its preference over other
explicit integration schemes. The reverse mode approach for accelerating the
computations applies immediately, but we do acknowledge that general second-
order sensitivity analysis for multi-stage RK quickly gets (notationally)
complicated. In part for the sake of time, but mostly to keep the paper self
contained, we have kept the Euler integration used in the results to be
consistent with the derivations.

While we had not previously considered implicit integration schemes, we do see
areas of improvement where we generalize our methods to them. In fact, the
authors in [ImplicitDDP] have nicely explained how to carry out DDP using
implicit dynamics, and there is clear opportunity to generalize this work
therein. However, those extensions are out of scope for this current paper.

We have added a few comments in the conclusions to point toward these
extensions, and these are highlighted in the next response item.

We also thank the reviewer for pointing out the Acado resource - this might
prove useful in our line of work.

n the plus side, and this is unfortunately not mentioned, their approach
should also work for other numerical optimal control methods such as multiple
shooting. It would be great to add a short discussion about this.

We agree with the reviewer that the contributions of our work are not only
specific to single-shooting DDP and that the methods are immediately
applicable to multi-shooting DDP variants. With this in mind, we revised the
‘LABEL:sect:_Conclusions’ as follows: [2.5em]0em  While AD tools are
convenient, they are general purpose, and thus may not be optimal. Alternative
analytical methods for taking derivatives of rigid-body dynamics can
accumulate the derivatives recursively, as in [carpentier2018analytical,
lee2005newton, sohl2001recursive, jain1993linearization]. Recently, we
extended the modified RNEA algorithm with an analytical accumulation of its
first-order partials in a reverse-mode fashion, and further evaluation of this
result is of immediate interest. Moreover, we aim to extend this work to
address rigid-body dynamics with contacts [DDPwithContacts, li2020hybrid] by
using similar approaches as in this paper. This generalization would then
allow for use with hybrid dynamic systems that arise in legged locomotion
problems.

There are many other opportunities that this work motivates as next steps.
While we noted that the presented work is general for any explicit integration
scheme, we also see opportunity to extend this work for implicit integration
and implicit DDP [ImplicitDDP]. Further, while the DDP used here was a single-
shooting solver, our contributions could be used in multi-shooting DDP and
other numerical optimal control solvers.

inor comments: - typo ”complied C/C$++$” $->$ ”compiled C/C$++$.”?

This typo has been removed. Thank you for bringing it to our attention.

verall the paper is decently written and is insightful. Its practical use is
possibly slightly lower than one might hope for initially (but maybe a good
example can be provided in the future). The work and analysis itself is
carried out thoroughly, making it a decent conference paper.

We thank the reviewer for their positive comments in regards to this work.

##  7 Reviewer 8

his paper provides a solid contribution, presented clearly and thoughtfully.
It is a nice and natural extension of the work on efficient differentiation of
manipulator equations. I am not fully convinced that it will have a big
impact, because it requires more effort to implement (especially the ABA
versions) than iLQR, and even if the convergence is quadratic vs super-linear,
it’s still not clear that DDP always converges better than iLQR. I think the
basin of attraction is probably more important than convergence rate near
optimality.

We thank the reviewer for their comments. To ease implementation effort and
for ease of understanding, we provide our Matlab code at
https://tinyurl.com/468ynkuu. We have plans to implement our derivatives in
Pinocchio as a next step to enable this work to achieve its full potential
impact on the community.

Moreover, towards making a connection to a basin of attraction, we randomize
the initial controller for the $7-$link KUKA LBR manipulator and present the
following findings along with Fig. 2

[2.5em]0em  Since the solution of an OCP is dependent on the initial
conditions, we randomize the initial controller using Ornstein-Uhlenbeck
process noise for the $7-$link KUKA LBR manipulator and solved the OCP problem
using either iLQR/DDP in $40$ separate instances. The first pane of Fig.
LABEL:fig:RandomCtrl illustrates the progression of that manipulator from an
initial configuration to the balanced upright configuration along an optimal
trajectory. Subplot (a) of Fig. LABEL:fig:RandomCtrl illustrates the empirical
probability density function (pdf) of the number of iterations over those
optimizations. Subplot (b) illustrates the pdf of the log of the final cost of
each optimal solution. As illustrated in subplot (a), in $95\%$ of the
simulations, iLQR had a higher number of iterations. On average, iLQR had
three times as many iterations as DDP. From subplot (b), we note that most of
the simulations regardless of DDP/iLQR converged to similar solution; in fact,
DDP converged to a different solution than iLQR in only three instances.
Figure LABEL:fig:RandomCtrl illustrates that the inclusion of second-order
information in DDP will result in similar converged solution and in fewer
iterations. This result is powerful in that the addition of second-order
information results in algorithm that converges in less iterations and to the
same minima.

Figure 2: Simulations on $7-$link KUKA LBR manipulator considering randomized
initial controller given by an Ornstein-Uhlenbeck process noise.

he major question / request I have for the paper is to clarify the role of the
manipulator equations in the contribution. First off, one can read the title
and the abstract, and have no idea that some of your main contributions are
specific to recursive multibody dynamics. Can you please make that clear?

To better capture our contributions and the domain specific application of our
work, we re-titled this paper as ‘Accelerating Second-Order Differential
Dynamic Programming for Rigid-Body Systems.’ We note that the first
contribution of using reverse-mode to calculate the necessary quantities for
DDP is a general contribution, and is not specific to rigid-body systems. But
the second contribution where we extend the relation of the first-order
sensitivity of $\operatorname{ID}$ to $\operatorname{FD}$ to the second-order
case is specific to rigid-body systems. We added a sentence to the abstract to
make this clear as follows: [2.5em]0em We also show how the structure of the
dynamics can be used to further accelerate these computations in rigid-body
systems.

econd, I think the initial contribution, in IIIA, where you multiply by the
lagrange multiplier before taking the second derivative, does *not* actually
seem specific to multibody systems. Unless I”m missing something, I would
recommend that you write that first contribution in the more general form,
e.g. d(lambda*fx)dx, instead of using FD already? That is quite a useful
observation for DDP more generally. I’m surprised nobody has used it
before(?), but I’ll certainly remember it now. Very nice. It’s quite
reminiscent of the way one computes Jdot*qdot jointly instead of ever
computing Jdot by itself when doing task-space derivatives. The numerical
evaluations are thoughtful and I think sufficient. It does leave one with the
impression that iLQR is probably just fine.

We thank the reviewer for their analysis and do note that indeed the first
contribution is general whereas the second contribution is specific to rigid-
body systems. With this in mind, we rewrite the contributions section as
follows: [2.5em]0em  This work presents and combines several advances for
reducing the computational complexity of computing second-order dynamics
sensitivity information in DDP. The final result is a method for computing
this information with the same computational complexity as first-order
dynamics derivatives (e.g., as in iLQR). The contributions are (I) the use of
reverse-mode automatic differentiation (AD) [grievank2000principles] to
compute second-order derivatives needed in DDP. This contribution is general
to discrete-time dynamic systems and enables computation reductions compared
to methods that explicitly evaluate a derivative tensor for the dynamics.
Further, we show (II) how second-order information related to the forward
dynamics can be related to associated information from the inverse dynamics,
akin to first-order results in [carpentier2018analytical,
jain1993linearization]; this contribution is specific to rigid-body models.
Lastly, we (III) introduce a modification to the Recursive-Newton-Euler
Algorithm (RNEA) that supports this process and further reduces computational
demands. Figure LABEL:fig:Fig1 overviews a benchmark of the proposed methods
against iLQR and against DDP approaches that explicitly compute derivative
tensors.

Regarding whether anyone has used our main AD tricks before, we are not aware
of any that have done so for DDP. However, the use of reverse mode is quite
common for Hessian calculations. That is to say that many AD packages give you
the option to use reverse mode (followed by forward mode) when computing the
Hessian of a generic scalar function. In this sense, the trick of accelerating
Hessian computations is not specific to functions of the form,
$\bm{\lambda}^{T}\mathbf{f}(\mathbf{x})$, but rather to any scalar function
$g(\mathbf{x})$ with $\mathbf{x}\in\mathbb{R}^{n}$. Using reverse mode on
$g(\mathbf{x})$ followed by forward mode on $\nabla g(\mathbf{x})$ will always
enable computing a Hessian with cost proportional to $n$ times that of
computing $g(\mathbf{x})$ itself.

While the applicability of using reverse mode for Hessian evaluations was
clear from reading AD toolbox manuals, we weren’t able to find a reference for
genesis of that idea. Still, we have acknowledged it in the revised manuscript
as follows: [2.5em]0em This general strategy also applies to any dynamic
system $\mathbf{x}_{k+1}=\mathbf{f}(\mathbf{x}_{k},\mathbf{u}_{k})$ for
computing
$\frac{\partial}{\partial\mathbf{x}}[\bm{\lambda}^{T}\mathbf{f}_{\mathbf{x}}]=\bm{\lambda}\cdot\mathbf{f}_{\mathbf{x}\mathbf{x}}$
in DDP. While the use of reverse-mode accumulation to compute Hessians is a
standard option in AD packages, its use to accelerate DDP here is new.

We genuinely thank the reviewer for helping us to better communicate the full
scope of our contribution.

##  9 Reviewer 9

he paper analyzes DDP methods when applied to robots, and more in general to
rigid-body systems. It shows that computation times can be reduced by
exploiting the structure of the Euler-Lagrange equations, and cheaply evaluate
the second-order derivatives in the DDP algorithm. This gives computation
times for DDP comparable with iLQR, which is a first order method. Overall, I
think the quality of the paper is good. The results are relevant for the
robotics community, the presentation is rigorous but not too heavy, and the
paper is well written. I only have two doubts that I’d like the authors to
comment on. One is quite high level. The second concerns more directly the
paper and is more important. 1\. I find it odd that in robotics the complexity
of algorithms is studied in terms of asymptotic dependence on the number of
DoFs. Personally, I’ve never seen a robot with more than 30 DoFs. I think that
most of the times the constant that multiplies the power of n is as important
as the power of n itself! In this sense, the results in the paper are
convincing in showing that the proposed methods have better asymptotic
complexity than traditional DDP, but then we see that traditional DDP
outperforms all the competitors when n $\leq$ 7 (Figs. 1 and 5). In practice,
I’d say that most of the robots do not have more than 7 DoFs, no?

We are mostly interested in legged robots and we aim to apply our work to
quadrupeds. In the whole-body MPC setting, legged robots are of high state
dimensions $(>=7)$. For example, [tassa2012synthesis] and [koenemann2015whole]
use DDP for a $22-$ and $25-$DoF humanoid respectively whereas [li2020hybrid]
use a $7-$DoF quadruped; In this setting, the quadruped is modeled using
simple models and is only modeled in 2D, it should noted in fact that a
complete model of the quadruped has 18 DoF [MiniCheetahThesis]. Towards that
aim, we consider the scaling of our work to be of paramount importance towards
future online whole-body optimal control strategies.

didn’t found the numeric results in Sec. 4 to be very convincing: - The
convergence criteria ($<10^{-9}$) seems very small, and it also strongly
favours DDP methods over iLQR, correct? Does this mean that if we were to set
the tolerance to a more reasonable value (say $10^{-5}$) iLQR would be much
faster than what you propose?

Even at the suggested convergence criteria, DDP would still have a fewer
number of iterations – this is illustrated in Fig. 3 of this response, whereby
we impose the suggested cutoff convergence criteria. Moreso, we include more
analysis on the KUKA LBR manipulator whereby we ran $40$ simulations with
randomized control inputs. These results are included in Fig.
LABEL:fig:RandomCtrl in the main paper whereby it is illustrated that iLQR at
most instances had twice as many iterations as DDP.

Figure 3: Suboptimality Decay vs Iterations for $7-$link KUKA LBR manipulator
with the suggested cut-off of $10^{-5}$ imposed.

itemize[]

Why is the modified RNEA always slower than RNEA in Fig. 5? Shouldn’t the
first be an improvement of the second?

”As shown, evaluation of the second-order partials by tensor contraction takes
the longest time.” This is not quite true for $n<=7?$ \- ”with the modified
RNEA outperforming RNEA for $n>=18$”. It seems to me that the purple and the
black lines cross at 21 and not at 18?

Overall, I’d like the authors to comment more on why algorithms that in theory
appear to be strictly better than others frequently perform worse

These two issues were due to subtle implementation error (that we had
previously missed) where extra overhead cost of the Modified RNEA was included
in the timing cost compared to RNEA calculation. Since we are using CasADi,
every interface between CaADi and MATLAB has some overhead associated with it,
and since Modified RNEA has more steps, we previously had multiple CasADi
function calls.

We have since corrected these errors so that all comparison are on equal
footing (just a single CasADi function call to obtain the derivatives for each
method). The revised paper now includes the reproduced timing figures and we
include them here as Fig. 4.

[Timing of the partials] [DDP/iLQR Timing]

Figure 4: Reproduced Timing Figures

don’t completely get why the number iterations in Fig. 7 are different when
you use different methods to compute the second-order derivatives. Shouldn’t
they be the same? And, if yes, why didn’t you just include two bars in the
plot (iLQR vs DDP)?

We agree with the reviewer that since DDP/iLQR variants solve the same cost
function staring from the same initial conditions and for the same model, the
number of iterations should be the same for each DDP variant, and the same for
each iLQR variant. In the original Fig. 7 (now deleted), we acknowledge that
the each DDP variant and each iLQR variant might have had slightly different
number of iterations. This difference was either one or two iterations for
either DDP/iLQR group. Due to machine precision, there may exist small
numerical differences that accumulate leading to minor differences in the
number of iterations. We have since deleted Fig. 6 and replicated the same
information in Fig. 9 We show the same result that iLQR variants take
significantly more iterations as compared to the DDP variants.

ere is a list of further minor comments:

  * •

Abstract: I would mention explicitly that the method applies to rigid-body
systems, and not to dynamical systems in general. Maybe just a word in the
first sentence. Maybe mention this also in Sec. 1A?

To clarity that our main contributions are specific to rigid-body systems, we
have retitled this paper as ‘Accelerating Second-Order Differential Dynamic
Programming for Rigid-Body Systems.’ We kindly note that our first
contribution of using reverse mode AD to accelerate DDP applies to dynamic
systems in general, and we have reworded our contributions in the introduction
to emphasize which contributions apply generally vs. only to rigid-body
systems. The ‘LABEL:sect:_Contributions’ section now reads as follows:

[2.5em]0em  This work presents and combines several advances for reducing the
computational complexity of computing second-order dynamics sensitivity
information in DDP. The final result is a method for computing this
information with the same computational complexity as first-order dynamics
derivatives (e.g., as in iLQR). The contributions are (I) the use of reverse-
mode automatic differentiation (AD) [grievank2000principles] to compute
second-order derivatives needed in DDP. This contribution is general to
discrete-time dynamic systems and enables computation reductions compared to
methods that explicitly evaluate a derivative tensor for the dynamics.
Further, we show (II) how second-order information related to the forward
dynamics can be related to associated information from the inverse dynamics,
akin to first-order results in [carpentier2018analytical,
jain1993linearization]; this contribution is specific to rigid-body models.
Lastly, we (III) introduce a modification to the Recursive-Newton-Euler
Algorithm (RNEA) that supports this process and further reduces computational
demands. Figure LABEL:fig:Fig1 overviews a benchmark of the proposed methods
against iLQR and against DDP approaches that explicitly compute derivative
tensors.

  * •

Sec. 1: Not sure that ”online optimization may allow a manipulator to rapidly
tailor its movements to a grasped object whose weight is not known a priori”
is a great example. I know what you mean, but it seems to imply that the
optimization is somewhat adaptive. I would rephrase.

We thank the reviewer for their comments regarding this sentence. For clarity,
we have updated this sentence to read as follows: [2.5em]0em For example,
online optimization may enable legged systems to tailor their gaits to sensed
terrains and inevitable disturbances or may enable manipulators to rapidly
synthesize efficient motions when transporting new objects.

  * •

Sec. 1: ”Recently, whole-body trajectory optimization is becoming …” I believe
it is worth mentioning Dai et al. ”Whole-body Motion Planning with Centroidal
Dynamics and Full Kinematics” here. This paper has been published before the
three papers you cite, and fits squarely in the description?

We agree with the reviewer’s sentiment and have added this citation in the
introduction.

  * •

Sec. 2B: I don’t think the ”Hamiltonian point of view” (Eq. (6) and the
equation below it) is self evident. I would cite something if possible, or
give a little more details.

We aim to clarify this by citing [mayne1966DDP] and noting that
$\bm{\lambda}\triangleq V_{\mathbf{x}}^{\prime}$. The revised statement reads
as follows: [2.5em]0em These coefficients (LABEL:eqn:coeffs) could
alternatively be viewed through derivatives of the Hamiltonian

$H_{k}(\mathbf{x},\mathbf{u},\bm{\lambda})=\ell_{k}(\mathbf{x},\mathbf{u})+\bm{\lambda}^{T}\mathbf{f}(\mathbf{x},\mathbf{u})\,,$
(1)

where $\bm{\lambda}\triangleq V_{\mathbf{x}}^{\prime}$ is the co-state vector.
For example,[mayne1966DDP]

$\displaystyle Q_{\mathbf{x}}$
$\displaystyle=\nabla_{\mathbf{x}}H_{k}(\mathbf{x}_{k},\mathbf{u}_{k},V_{\mathbf{x}}^{\prime})\quad{\textrm{
and}}$ $\displaystyle Q_{\mathbf{u}\mathbf{x}}$
$\displaystyle=\nabla^{2}_{\mathbf{u}\mathbf{x}}H_{k}(\mathbf{x}_{k},\mathbf{u}_{k},V_{\mathbf{x}}^{\prime})+\mathbf{f}_{\mathbf{u}}^{T}V_{\mathbf{x}\mathbf{x}}^{\prime}\mathbf{f}_{\mathbf{x}}\,.$

  * •

Sec. 2C: Please provide a citation when you say that DDP converges
quadratically and iLQR super linearly.

In [tassa2012synthesis], they note that ”iLQG is a variant of the classic
Differential Dynamic Programming (DDP) algorithm, the main difference being
that only first rather than second derivatives of the dynamics are used. This
means that iLQG no longer exhibits the quadratic convergence properties of
DDP.” In Section IIC we have added this reference as well.

Extremely minor comments:

  * •

Sec. 2A: Strictly speaking, I believe that calling C(q, qdot) the ”Coriolis
terms” is incorrect. I’d call them centrifugal and Coriolis terms or quadratic
velocity terms. (The second is taken from Siciliano et al. ”Robotics:
Modelling, Planning and Control” pages 256-257.)

We updated the language from “Coriolis terms” to “Coriolis and centrifugal
terms”.

  * •

\- Sec. 3: I think your contribution is very clear at this point, I wouldn’t
say ”Contribution” in the section title.

Thank you for the feedback. We updated the section title of Section III from
“Contribution: Efficient Computation of Second-Order Derivatives for DDP” to
“LABEL:sect:Derivs”.

  * •

Sec. 3A: “suppose it can be computed with complexity…” I’d say evaluated
instead of computed? Also, I don’t get the reason of the sentence ”where
C(.,.) is a function of its arguments.”

For clarity, we have edited this paragraph to read as follows: [2.5em]0em To
motivate our approach for accelerating the evaluation of (6), consider any
vector-valued function
$\mathbf{g}(\mathbf{x}):\mathbb{R}^{n}\rightarrow\mathbb{R}^{m}$. Given any
fixed vector $\bm{\gamma}\in\mathbb{R}^{m}$, reverse-mode derivative
accumulation [grievank2000principles] provides an approach to evaluate the
vector-Jacobian product

$\bm{\gamma}^{T}\frac{\partial\mathbf{g}}{\partial\mathbf{x}}=\sum_{j=1}^{m}\gamma_{j}\frac{\partial
g_{j}}{\partial\mathbf{x}}=\frac{\partial}{\partial\mathbf{x}}\left[\bm{\gamma}^{T}\mathbf{g}\right]$
(2)

with the same complexity as evaluating $\mathbf{g}(\mathbf{x})$. It is
implemented in many AD packages (e.g., CasADi [andersson2019casadi]). In
practice, the cost of evaluating (3) can be bounded above by a small constant
times the cost of evaluating $\mathbf{g}(\mathbf{x})$, and the constant is
often three to four [grievank2000principles, Section 3.3].

##  10 Reviewer 11

his paper presents a novel methodology for the optimization of trajectories in
multi-link robots. A modified DDP algorithm is presented to avoid the
computation of tensor-related variables which increase the computational
demand when a large number of links is considered. The major contribution are
very clear since the Introduction Section of the paper. Final simulation
results asses the performances of the proposed method. Moreover, the proposed
approach has been tested under stressing conditions (e.g., with dissipative
controllers as starting controllers) so that its performances are clearly
demonstrated. Organization and Style: The paper is well organized and its
style is adequate for the IEEE Robotics and Automation Letters journal.
Technical Accuracy: The proposed methodology is presented with a high level of
accuracy. Mathematical details are rigorous and the technical arguments are
sound and solid. Presentation: The paper is well presented. Just few remarks:
It would be nice to see a comparison between the actual trajectories obtained
with the traditional iLQR and the ones obtained via the tensor-free DDP
approach.

Due to the page constraints, we chose to omit the trajectories in the main
paper but include them (see Fig. 5) here for your convenience. This figure
encodes much of the same information as the snapshots in the first panel of
Fig. 6.

Figure 5: Trajectories produced by DDP for the KUKA LBR manipulator. Subplot
(a) shows the state trajectory whereby $q_{1}$ settles to $\frac{\pi}{2}$
since $\mathbf{q}_{1}$ is measured relative to the world horizontal axis. The
other states $q$ and $\dot{q}$ settle to zero since each state is measured
relative to the previous joint. Subplot (c) shows the torque of each joint and
finally subplot (d) shows the sub-optimality decay

We interpret that at the heart of the reviewer’s question is concern regarding
the quality of the trajectories obtained by DDP/iLQR. Towards that aim, we
randomize the initial controller for the $7-$link KUKA LBR manipulator and
present the following findings along with Fig. 6 [2.5em]0em  Since the
solution of an OCP is dependent on the initial conditions, we randomize the
initial controller using Ornstein-Uhlenbeck process noise for the $7-$link
KUKA LBR manipulator and solved the OCP problem using either iLQR/DDP in $40$
separate instances. The first pane of Fig. LABEL:fig:RandomCtrl illustrates
the progression of that manipulator from an initial configuration to the
balanced upright configuration along an optimal trajectory. Subplot (a) of
Fig. LABEL:fig:RandomCtrl illustrates the empirical probability density
function (pdf) of the number of iterations over those optimizations. Subplot
(b) illustrates the pdf of the log of the final cost of each optimal solution.
As illustrated in subplot (a), in $95\%$ of the simulations, iLQR had a higher
number of iterations. On average, iLQR had three times as many iterations as
DDP. From subplot (b), we note that most of the simulations regardless of
DDP/iLQR converged to similar solution; in fact, DDP converged to a different
solution than iLQR in only three instances. Figure LABEL:fig:RandomCtrl
illustrates that the inclusion of second-order information in DDP will result
in similar converged solution and in fewer iterations. This result is powerful
in that the addition of second-order information results in algorithm that
converges in less iterations and to the same minima.

Figure 6: Simulations on $7-$link KUKA LBR manipulator considering randomized
initial controller given by an ornstein-Uhlenbeck process noise.

would spend a little bit more words on the DDP theory in general.

Due to the page limitation, we only chose to focus on the most salient
features of DDP/iLQR. While our discussion of DDP is short, we refer the
reader to a multitude of DDP/iLQR literature [mayne1966DDP, li2004iterative,
tassa2012synthesis, liao1991convergence, koenemann2015whole].

would introduce the “AD” acronym before citing Figure 2, I think now it is a
bit confusing.

In this revised paper, we ensure that the term ‘Automatic Differentiation
(AD)’ is introduced prior to the reference of Fig. 2.

We would appreciate feedback from the reviewer as to whether the reference to
Figures 2 and 3 appeared too early, or if they were helpful as an early
roadmap to the development.

s the iLQR used in the comparison optimized? Some parts of the iLQR are
optimized for the computation? Please provide more info on this aspect.

Given a cost function, both iLQR and DDP optimize the cost function while
still respecting the dynamics constraints. Following the initial review, DDP
and iLQR were built on the same code base with the only difference between
them is that DDP accounts for the the tensorial terms (e.g.,
$\mathbf{f}_{\mathbf{x}\mathbf{x}}$ in different manners) and adds
regularization if $Q_{\mathbf{u}\mathbf{u}}$ is not positive definite at any
point. Since DDP and iLQR are built on the same code base, in implementation,
the use of DDP or iLQR is controlled by a binary flag that can be set in the
code. This revised design represents our best-faith effort at fairness in the
comparison.

ow your method behaves in the neighbors of kinematic singularities? Please
provide more info on this aspect.

We humbly acknowledge that due to the limited revision time of RA-L and the
numerous reviews to this document, we did not have the time to simulate a
manipulator with a kinematic singularity as its end goal. We do wish to
highlight, however, that we have no Jacobian inversions (or pseudoinversions
for that matter) in our implementation, such that the code should be
numerically robust to these cases.

t would be very useful (but it is not mandatory for the acceptance of the
paper, just an advice) if you share some code so that such method can be
easily implemented an replicated for other applications and works.

The MAtlab implementation for this work has been added online to github at:
https://tinyurl.com/468ynkuu. This is an early release and we will look to
improve it before any final submission.

##  12 Reviewer 12

he paper proposes the use of ”reverse-mode derivative accumulation” [27] to
reduce the computational complexity of second-order dynamics sensitivity
information in Differential Dynamic Programming (DDP), from
$\mathcal{O}(n^{3})$ to $\mathcal{O}(n^{2})$, matching that of iLQR, by
eliminating the expensive tensor product. To further accelerate the
computation, it makes use of the duality between forward dynamics (FD) and
inverse dynamics (ID), and uses RNEA instead of ABA. While general, these
improvements appear to be motivated by the intention to obtain high-fidelity
trajectories for floating-base, high-DoF systems–specifically quadrupeds.
Overall, the paper is well-written. The sections build upon each other to form
a coherent story. The introduction covers the most-relevant DDP-based
literature. An extensive amount of mathematics is thoroughly discussed and
appears to be correct. The following comments can be made: 1) The lack of
”Conclusion” section is apparent. A reader who is unfamiliar with the topic
will find the material challenging to understand. Their understanding should
be reinforced at the end of the text.

We have updated the last section such that the first paragraph serves as the
conclusion and the remaining paragraphs serve as the ‘future directions’ of
this work. With this refocus of the section, we reinforce the presented work
at end of the paper and ensure that future directions of this work are clear.
The last section is structured as follows: [2.5em]0em  This work has made use
of reverse-mode AD tools to quickly evaluate different approaches for
obtaining the derivatives needed by DDP. We also extended the relationship
between the first-order sensitivities of $\operatorname{ID}$ and
$\operatorname{FD}$ [carpentier2018analytical, jain1993linearization] for
rigid-body systems to the second-order case. The combination of this new
approach and AD tools allows for the evaluation of the needed derivatives in
DDP with the same complexity as iLQR. Lastly, we introduce a restructuring of
RNEA to derive a modified RNEA that returns $\mu^{T}\operatorname{ID}$ in
$\mathcal{O}(n)$ complexity, and enables the fastest DDP algorithm.

he introduction could explicitly reference other direct trajectory
optimization methods like collocation, transcription, or multiple shooting
(e.g.,

  * •

$[*],[^{]}).[*]$ M. Diehl, H. G. Bock, H. Diedam, and P.B. Wieber, ‘Fast
direct multiple shooting algorithms for optimal robot control’, in Fast
motions in biomechanics and robotics: Optimization and feedback control, M.
Diehl and K. Mombaur, Eds. Berlin, Heidelberg: Springer Berlin Heidelberg,
$2006,pp.65–93.[^{]}$

  * •

M. Kelly, ‘An Introduction to Trajectory Optimization: How to Do Your Own
Direct Collocation’, SIAM Rev., vol.
$59,no.4,pp.849–904,Jan.2017,doi:10.1137/16M1062569.$

We thank the reviewer for providing additional citations that strengthen this
paper. We have since added the recommended citations to the introduction.

he contribution section should cite [grievank2000principles] and could make
use of bullet points.

The ‘LABEL:sect:_Contributions’ section now cites [grievank2000principles],
uses enumeration, and reads as follows: [2.5em]0em  This work presents and
combines several advances for reducing the computational complexity of
computing second-order dynamics sensitivity information in DDP. The final
result is a method for computing this information with the same computational
complexity as first-order dynamics derivatives (e.g., as in iLQR). The
contributions are (I) the use of reverse-mode automatic differentiation (AD)
[grievank2000principles] to compute second-order derivatives needed in DDP.
This contribution is general to discrete-time dynamic systems and enables
computation reductions compared to methods that explicitly evaluate a
derivative tensor for the dynamics. Further, we show (II) how second-order
information related to the forward dynamics can be related to associated
information from the inverse dynamics, akin to first-order results in
[carpentier2018analytical, jain1993linearization]; this contribution is
specific to rigid-body models. Lastly, we (III) introduce a modification to
the Recursive-Newton-Euler Algorithm (RNEA) that supports this process and
further reduces computational demands. Figure LABEL:fig:Fig1 overviews a
benchmark of the proposed methods against iLQR and against DDP approaches that
explicitly compute derivative tensors.

ec. II is titled ”Background: Trajectory Optimization” but only mentions DDP.

We agree with the reviewer that we only focus on DDP/iLQR within this paper
and that there exists numerous trajectory optimization frameworks such as
collocation, multi-shooting etc. For clarity, we retitled Section II as
‘LABEL:sect:_DDPBackground.’

ertain subsections (e.g., III-A, III-B) excessively introduce the subsequent
subsection. The authors should omit such paragraphs or place them inside the
respective sections so that they can be read independently.

We have cut back on introducing the subsequent subsection where it might not
be necessary to do so.

ll figures have too much negative vspace between the graph and the captions
which negatively impacts reading.

In light of this of review, we have increased the amount of white space around
all figures. We hope that the spacing is adequate and enables for clear
reading.

here is no need to reference Fig. 2 in Sec. III. It diverts the reader’s
attention before they can comprehend significant of the different approaches
and understand the terms and acronyms in the diagram.

We have removed the early references to these figures and placed them at the
end of their respective subsections. For example, Fig. 2 is now referenced
near the end of III.A and Fig. 3 is referenced after the main derivation in
III.B. By delaying these references, we believe that this will minimize
confusion, while still enabling the figures to reinforce the main take-aways
of the section.

We acknowledge that the reader does not have all information evaluate 100% of
Figure 2 by the end of Section III.A, but we want to avoid the reader having
multiple pages straight of mathematics to follow without any visual aids.
Recognizing that this could be a matter of personal taste, we have requested
feedback from all reviewers on this point (at the opening of this document)
and will look to follow the suggestions of the majority for any future
versions of this manuscript.

ec. III-A introduces ”reverse-mode accumulation” which boils down to a
function and its gradient having the same computational complexity; this is
said in a long-winded/inverted way.

We clarify our introduction of reverse-mode accumulation as follows below. We
would welcome any further feedback on whether we should strive to make it more
concise. [2.5em]0em To motivate our approach for accelerating the evaluation
of (6), consider any vector-valued function
$\mathbf{g}(\mathbf{x}):\mathbb{R}^{n}\rightarrow\mathbb{R}^{m}$. Given any
fixed vector $\bm{\gamma}\in\mathbb{R}^{m}$, reverse-mode derivative
accumulation [grievank2000principles] provides an approach to evaluate the
vector-Jacobian product

$\bm{\gamma}^{T}\frac{\partial\mathbf{g}}{\partial\mathbf{x}}=\sum_{j=1}^{m}\gamma_{j}\frac{\partial
g_{j}}{\partial\mathbf{x}}=\frac{\partial}{\partial\mathbf{x}}\left[\bm{\gamma}^{T}\mathbf{g}\right]$
(3)

with the same complexity as evaluating $\mathbf{g}(\mathbf{x})$. It is
implemented in many AD packages (e.g., CasADi [andersson2019casadi]). In
practice, the cost of evaluating (3) can be bounded above by a small constant
times the cost of evaluating $\mathbf{g}(\mathbf{x})$, and the constant is
often three to four [grievank2000principles, Section 3.3].

ater, also in Sec. III-A, in the sentence ”Since FD can be calculated (…)” it
says that $\eta^{T}$… can be computed in $\mathcal{O}(n)$ but in Fig. 2 it
appears to be incorrectly stated as $\mathcal{O}(n^{2})$.

We thank the reviewer for noting this error in our graphic. We have corrected
the figure and ensured that the rest of the graphics are correct as well. We
include the corrected figure below Fig. 7 for your convenience.

Figure 7: Computational approach for the partials needed in _iLQR via ABA_ ,
_DDP via ABA_ , and _DDP via Tensor Contraction_ methods

he sentence ”(…) simple relationship [21], [22].” should end in a semicolon or
the paragraph should be rephrased.

This sentence has been updated to read as follows: [2.5em]0em  Derivatives of
$\operatorname{FD}$ are, however, analytically related to the derivatives of
$\operatorname{ID}$ by [carpentier2018analytical, jain1993linearization]

$\displaystyle\left.\frac{\partial\operatorname{FD}}{\partial\mathbf{z}}\right|_{\mathbf{q}_{0},{\dot{\mathbf{q}}}_{0},\bm{\tau}_{0}}$
$\displaystyle=-\mathbf{M}^{-1}(\mathbf{q}_{0})\left.\frac{\partial\operatorname{ID}}{\partial\mathbf{z}}\right|_{\mathbf{q}_{0},{\dot{\mathbf{q}}}_{0},\bm{\tau}_{0}}\quad{\textrm{and}}$
(4)
$\displaystyle\left.\frac{\partial\operatorname{FD}}{\partial\bm{\tau}}\right|_{\mathbf{q}_{0},{\dot{\mathbf{q}}}_{0},\bm{\tau}_{0}}$
$\displaystyle=\mathbf{M}^{-1}(\mathbf{q}_{0})=\left(\left.\frac{\partial\operatorname{ID}}{\partial{\ddot{\mathbf{q}}}}\right|_{\mathbf{q}_{0},{\dot{\mathbf{q}}}_{0},{\ddot{\mathbf{q}}}_{0}}\right)^{-1}\,,$
(5)

where $\mathbf{z}$ can either represent the vector $\mathbf{q}$ or
${\dot{\mathbf{q}}}$, and where
${\mathbf{q}_{0},{\dot{\mathbf{q}}}_{0},{\ddot{\mathbf{q}}}_{0}}$, and
$\bm{\tau}_{0}$ refer to the point of linearization.

itemize

In Sec. III-C it is stated ”(…) the new approach gives lower computational
cost since RNEA is simpler than ABA.”; however, this is not reflected by the
first-order partials in Fig. 4. A Minor benefit can be observed in the case of
second-order partials for DoF $>6$. The authors should explain why this is the
case.

In Sec. IV-A, the statement ”Above n = 6, the second-order partials via
RNEA/modified RNEA were faster than second-order partials via ABA (…)” is also
incorrect. In Fig. 4, at n = 6, only RNEA, and not Modified RNEA, is faster
than ABA. This should be corrected.

In Sec. $IV-A$, the statement ”$(...)1.5$ times more computation time than
first order partials for $n\geq 10$” is incorrect. In Fig. $4$, the ratio at
$n=10$ is closer to $3.3$. This should be corrected.

This issue was due to subtle implementation error (that we had previously
missed) where extra overhead cost of the Modified RNEA was included in the
timing cost compared to RNEA calculation. Since we are using CasADi, every
interface between CaADi has some overhead associated with it, and since
Modified RNEA has more steps, we previously had multiple CasADi function
calls.

We have since corrected these errors so that all comparison are on equal
footing (just a single CasADi function call to obtain the derivatives for each
method). We include the reproduced timing figures as Fig. LABEL:fig:Derivs and
Fig. LABEL:fig:LogTime in the main paper. We attach those figures here as Fig.
8 for your convenience.

[Timing of the partials] [DDP/iLQR Timing]

Figure 8: Reproduced Timing Figures

ig. 5 shows all methods taking a substantial amount of time irrespective of
the DoF. This indicates that the implementation (MATLAB) is inefficient. It
raises the question: To what degree has the implementation affected the
comparison and would the ratios be the same in a more efficient
implementation?

Following the initial review, DDP and iLQR were built on the same code base
with the only difference between them is that DDP accounts for the the
tensorial terms (e.g., $\mathbf{f}_{\mathbf{x}\mathbf{x}}$ in different
manners) and adds regularization if $Q_{\mathbf{u}\mathbf{u}}$ is not positive
definite at any point. Since DDP and iLQR are built on the same code base, in
implementation, the use of DDP or iLQR is controlled by a flag that can be set
in the code. This revised design represents our best-faith effort at fairness
in the comparison.

Currently, this work is implemented in Matlab but the $\mathcal{O}(\cdot)$
complexity quoted for the methods is irrespective of the implementation
platform. With this in mind, the value of the work is in its relative
comparison between the methods. The absolute timing will further be platform
dependent. Nonetheless, we eventually plan to implement this work in C++ with
our new derivatives computed via extensions to Pinocchio. Based on some of our
recent work, we would expect the derivative calculation times to be multiple
orders of magnitude faster in C++, due to a combination of efficient code, and
special purpose dynamics derivatives routines, as opposed to general purpose
auto-differentiation tools. However, we stand by the value of our current
work, since the flexibility of AD allowed us to test many configurations and
to do so in a fair manner.

he computation times shown in Fig. 6 are extremely high, regardless of the
fact that the robot has 15 links. Considering that the paper is framed around
the idea of online trajectory optimization, this is not a persuasive result.
Even if the proposed approach accelerates DDP, it appears to be extremely slow
considering that other labs achieve 50-100 Hz legged locomotions using iLQR.

This work focuses on the introduction of our theoretical contribution to the
wider community. We too are motivated by the excellent work in other labs that
has led to the 50-100Hz MPC implementations. However, those cases are 1) in
C++ 2) often warm-started and 3) usually only run for a single iteration in
the MPC loop. While we humbly acknowledge that we have a ways to go on
practical implementation, since this work is implemented in Matlab, we hope
the reader will consider the results compared relative to each other and not
in absolute timing. We have added a note that the results are in Matlab to the
figure captions of key results to assist with their proper interpretation.

In the future, we do plan to implement this work in compiled C$++$ and in
hardware whereby we will be able to report more persuasive absolute timing
computations. In the interim, we have released source code of our AD-based
comparison on git at https://tinyurl.com/468ynkuu as well. This code is a
first release, and we look to improve it for usability before any final
publication.

dditional minor comments are:

  * •

”This approach in diagrammed” should be ”is”

Fixed!

  * •

In Sec. III-A, ”reverse mode” should by hyphenated.

We have since corrected any instance where ‘reverse-mode’ was not hyphenated.

  * •

In Sec. III-B, ”state of the art” should by hyphenated.

This term is now hyphenated.

  * •

Sec. IV-C, the sentence ”Succinctly, ($\ldots$)” adds no new information,
merely repeats the argument made in the previous sentence. Non-native speakers
are likely to not know the meaning of ”succinctly”.

This sentence has now been removed in the revised paper.

  * •

In the caption of Fig. 7, ”($\ldots$) for an 15-link ($\ldots$)” should be
”a”.

The caption has now been updated.

  * •

Fig. 6 and Fig. 7 should be merged into a single figure as subfigures.

In the revised paper, we deleted Fig. 7 and reintroduce the same information
as Fig. 7 as a subplot in the newly reproduced Fig. LABEL:fig:RandomCtrl.

  * •

Fig. 8 is unnecessarily large.

We have reduced the size of the figure.

  * •

name of the physics simulator should stated in the caption of Fig. 9

The caption of this figure (now Fig. 8 in the text) was updated to mention
that we used Matlab’s Robotic Toolbox to run the simulation.

  * •

No video has been included.

In lieu of a video, we have worked to provide an open implementation of our
work, which is available at https://tinyurl.com/468ynkuu. This is a first
release, and we will look to improve its usability before the final
submission.

##  13 Reviewer 13

his paper presents one main improvement and two secondary improvements to
reduce the computational cost of Differential Dynamic Programming (DDP)
algorithm. The main contribution is the use of a reverse-mode derivative
accumulation technique to evaluate DDP (usually with complexity $O(n^{3})$, n
the number of degrees of freedom) with the same computational complexity than
iterative Linear Quadratic Regulator (iLQR) $(O(n^{2}))$. It allows to have a
$O(n^{2})$ DDP and so to keep second-order dynamics information without
increasing complexity. Other contributions include a way to use Recursive
Newton Euler Algorithm (RNEA) to perform computations (simpler than
Articulated-Body Algorithm ABA) and a slight modification of the RNEA
algorithm to further reduce computational demands. The content is well
presented and is explicit enough to be easily followed. Contributions are
clearly outlined which is a good point. The bibliography is extensive and
supports the contributions.

Thank you for the comments regarding our work. We have taken several steps to
improve our work as follows:

ntroduction: I would rephrase the first sentence of ”Specific contributions”.
Second-order dynamics are already part of DDP, it’s more that your work
combines several advances to reduce its computational complexity to get
roughly the same than iLQR.

We thank the reviewer for noting the ambiguity in the sentence. The specific
contribution section now reads as follows: [2.5em]0em  This work presents and
combines several advances for reducing the computational complexity of
computing second-order dynamics sensitivity information in DDP. The final
result is a method for computing this information with the same computational
complexity as first-order dynamics derivatives (e.g., as in iLQR). The
contributions are (I) the use of reverse-mode automatic differentiation (AD)
[grievank2000principles] to compute second-order derivatives needed in DDP.
This contribution is general to discrete-time dynamic systems and enables
computation reductions compared to methods that explicitly evaluate a
derivative tensor for the dynamics. Further, we show (II) how second-order
information related to the forward dynamics can be related to associated
information from the inverse dynamics, akin to first-order results in
[carpentier2018analytical, jain1993linearization]; this contribution is
specific to rigid-body models. Lastly, we (III) introduce a modification to
the Recursive-Newton-Euler Algorithm (RNEA) that supports this process and
further reduces computational demands. Figure LABEL:fig:Fig1 overviews a
benchmark of the proposed methods against iLQR and against DDP approaches that
explicitly compute derivative tensors.

rajectory Optimization: In(3), V*(xk) “can be construed as”

Thank you for catching this typo. We have updated the sentence to “the value
function can be constructed as…”.

fficient computation: Define meaning of $\bm{\xi}$ and $\bm{\eta}$ in
$\bm{\lambda}$

We have revised this paper to include the definition of these quantities as
follows: [2.5em]0em This section considers how dynamics partials enter the
partials of the Hamiltonian (1). The co-state vector is partitioned as
$\bm{\lambda}=[\bm{\xi}^{T},\bm{\eta}^{T}]^{T}$ where $\bm{\xi}$ and
$\bm{\eta}$ are the co-states associated with $\mathbf{q}$ and
${\dot{\mathbf{q}}}$ respectively. Via (LABEL:eqn:DiscreteDyn), the
Hamiltonian is then:

$\\!\\!H_{k}(\mathbf{x},\mathbf{u},\bm{\lambda})=\ell_{k}(\mathbf{x},\mathbf{u})+\begin{bmatrix}\bm{\xi}\\\
\bm{\eta}\end{bmatrix}^{T}\\!\left(\mathbf{x}+h\begin{bmatrix}{\dot{\mathbf{q}}}\\\
\operatorname{FD}(\mathbf{q},{\dot{\mathbf{q}}},\bm{\tau})\end{bmatrix}\right).$

fficient computation: In (9), brief explanation of why
$\bm{\xi}^{T}\frac{\partial}{\partial\mathbf{q}_{i}}[\frac{\partial}{\partial\mathbf{q}_{j}}{\dot{\mathbf{q}}}]$
is equal to 0

This follows since ${\dot{\mathbf{q}}}$ is not a function of $\mathbf{q}$. and
thus its derivative is zero. We have added a clarification as follows:
[2.5em]0em the partials can be written as

$\displaystyle\frac{\partial^{2}H}{\partial q_{i}\partial q_{j}}$
$\displaystyle=\frac{\partial^{2}\ell}{\partial q_{i}\partial
q_{j}}+h\begin{bmatrix}\bm{\xi}\\\
\bm{\eta}\end{bmatrix}^{T}\frac{\partial}{\partial
q_{i}}\left(\frac{\partial}{\partial
q_{j}}\begin{bmatrix}{\dot{\mathbf{q}}}\\\
\operatorname{FD}\end{bmatrix}\right)$
$\displaystyle=\frac{\partial^{2}\ell}{\partial q_{i}\partial
q_{j}}+h\bm{\eta}^{T}\frac{\partial}{\partial
q_{i}}\left[\frac{\partial}{\partial q_{j}}\operatorname{FD}\right]\,,$ (6)

where the simplification occurs since $\partial{\dot{\mathbf{q}}}/\partial
q_{j}=0$.

fficient computation:In Fig 2 the result of Reverse AD (top middle) should be
$O(n)$ and not $O(n^{2})$ since the point of reverse AD is not increasing the
complexity

We thank the reviewer for noting this error in the graphic. We have since
fixed the error and include here (see Fig. 9) for your convenience.

Figure 9: Computational approach for the partials needed in _iLQR via ABA_ ,
_DDP via ABA_ , and _DDP via Tensor Contraction_ methods

fficient computation: Clarify the sentences before ”We denote this method as
iLQR via RNEA.” since they mention ”calls to ABA” without explicitly bringing
up RNEA (it is used to get ID) State clearly that indirectly applying the mass
matrix inverse via calls to ABA allows to remain at $\mathcal{O}(n^{2})$
otherwise it is not clear why it is $\mathcal{O}(n^{2})$ and not
$\mathcal{O}(n^{3})$ in Fig 3 ”iLQR via RNEA”. Since as you pointed out other
state of the art methods have $\mathcal{O}(n^{3})$ to compute that.

The following paragraph now reads as follows: [2.5em]0em  These relationships
provide alternate methods for the first-order partials of $\operatorname{FD}$,
shown on the left side of Fig. LABEL:fig:DerivsRNEAgraph. The term
$\frac{\partial\operatorname{ID}}{\partial\mathbf{z}}$ can be computed in
$\mathcal{O}(n^{2})$ complexity with AD tools or specialized algorithms (e.g.,
[carpentier2018analytical, lee2005newton, sohl2001recursive,
jain1993linearization]). The explicit computation of the mass matrix inverse
can be avoided in (4) by instead applying it indirectly via $n$ calls to the
ABA algorithm ($\mathcal{O}(n)$) with the columns of
$\frac{\partial\operatorname{ID}}{\partial\mathbf{z}}$ as inputs for
$\bm{\tau}$. This approach evaluates (4) with total complexity
$\mathcal{O}(n^{2})$. Since (4) relies on RNEA to obtain $\operatorname{ID}$,
we denote this method as _iLQR via RNEA_. As an alternate approach, the
partials of $\operatorname{FD}$ are computed via (4) with $O(n^{3})$
complexity as follows. The partials of $\operatorname{ID}$ are computed with
$O(n^{2})$ complexity, the mass matrix inverse is computed once with
$O(n^{2})$ complexity (e.g., via [carpentier2018Inverse]), and a dense matrix-
matrix multiply in (4) finally sets the complexity at $O(n^{3})$. Since matrix
multiplications are optimized on modern hardware, this approach
[carpentier2018analytical] can be faster than the lower-order one. We next
extend the relation (4) to the second-order case.

fficient computation: Explain or add a reference to make it more clear why you
state that $M^{-1}\bm{\eta}$ is constant or it does not depend on q anymore
(since M(q)) At the end of III.C.

This point is subtle and one we have struggled with in our own communications
together. We have done our best to make it more clear as follows. [2.5em]0em
Once $\mathbf{q}$ and $\bm{\eta}$ are given, we fix the value of
$\bm{\mu}\triangleq\mathbf{M}^{-1}\bm{\eta}=\operatorname{FD}(\mathbf{q},0,\bm{\eta},0)$
(i.e., once fixed it no longer depends on $\mathbf{q}$).  [2.5em]0em Note that
the key reorganization $\frac{\partial}{\partial
q_{j}}[\bm{\mu}^{T}\frac{\partial\operatorname{ID}}{\partial
q_{i}}]=\bm{\mu}^{T}\frac{\partial^{2}\operatorname{ID}}{\partial
q_{i}\partial q_{j}}$ is possible since $\bm{\mu}$ is treated as constant
after its initial evaluation.

To be more notationally verbose about it we have that:

$\bm{\eta}^{T}\mathbf{M}^{-1}\frac{\partial^{2}\operatorname{ID}}{\partial
q_{i}\partial
q_{j}}=\left.\bm{\mu}^{T}\frac{\partial^{2}\operatorname{ID}}{\partial
q_{i}\partial
q_{j}}\right|_{\bm{\mu}=\mathbf{M}(\mathbf{q})^{-1}\bm{\eta}}=\left.\left(\frac{\partial}{\partial
q_{i}}\left[\bm{\mu}^{T}\frac{\partial\operatorname{ID}}{\partial
q_{j}}\right]\right)\right|_{\bm{\mu}=\mathbf{M}(\mathbf{q})^{-1}\bm{\eta}}$

So, the main idea is that we “fix” $\bm{\mu}$ as not depending on $\mathbf{q}$
during the process of taking derivatives. While our presentation approach in
the paper may not be as precise as above, noting the main idea of fixing
$\bm{\mu}$ avoids the lengthier equations as above and streamlines their
presentation. Since the equations in this section already don’t fit on one
line, we thought it was important to keep them simple as much as possible. We
hope this helps to clarify our choices, and we would certainly welcome and
value any further feedback from the reviewer on improving clarity for this
tricky point.

fficient computation: For the other derivatives, I do not understand where the
$–$ signs in front of $\Psi^{T}$ and $\Xi^{T}$ come from since
$\frac{\partial}{\partial q}(\eta^{T}\Psi)=\frac{\partial}{\partial
q}(\Psi^{T}\eta)=\Psi^{T}\frac{\partial\eta}{\partial
q}+[…]=\Psi^{T}\frac{\partial}{\partial q}(M\mu)-[\ldots]$. Same for $\Xi^{T}$

The partial of a matrix inverse is given as follows for any invertible matrix
$\mathbf{M}$ and any $q_{j}$ as

$\frac{\partial}{\partial
q_{j}}\mathbf{M}^{-1}=-\mathbf{M}^{-1}\frac{\partial\mathbf{M}}{\partial
q_{j}}\mathbf{M}^{-1}.$

We noted where this property in used in the current derivation: [2.5em]0em
where we used that $\partial\mathbf{M}^{-1}/\partial
q_{j}=-\mathbf{M}^{-1}(\partial\mathbf{M}/\partial q_{j})\mathbf{M}^{-1}$

We will repeat the derivation here for one of the other partials to show that
the minus sign appears similarly there. In terms of explaining what is going
wrong with the mathematics you presented, please note that $\bm{\eta}$ is the
co-state, which doesn’t depend on $\mathbf{q}$. The derivation remains tricky
since
$\frac{\partial}{\partial\mathbf{q}}\left(\bm{\eta}^{T}\bm{\Psi}\right)=\bm{\eta}\cdot\frac{\partial\bm{\Psi}}{\partial\mathbf{q}}$
where $\bm{\Psi}$ is a matrix, so
$\frac{\partial\bm{\Psi}}{\partial\mathbf{q}}$ is a tensor. So, to carry out
the derivative, we proceed element by element, and then will re-express the
final result back together in matrix form at the end.

Since
$\frac{\partial\operatorname{FD}}{\partial\dot{q}_{i}}=-\mathbf{M}(\mathbf{q})^{-1}\frac{\partial\operatorname{ID}}{\partial\dot{q}_{i}}$,
we have:

$\displaystyle\frac{\partial}{\partial
q_{j}}\left[\bm{\eta}^{T}\frac{\partial\operatorname{FD}}{\partial\dot{q}_{i}}\right]$
$\displaystyle=\frac{\partial}{\partial
q_{j}}\left[-\bm{\eta}^{T}\mathbf{M}^{-1}(\mathbf{q})\left.\frac{\partial\operatorname{ID}}{\partial\dot{q}_{i}}\right|_{\mathbf{q},{\dot{\mathbf{q}}},{\ddot{\mathbf{q}}}=\operatorname{FD}(\mathbf{q},{\dot{\mathbf{q}}},\bm{\tau})}\right]$
(7) (Using Product Rule and Chain Rule…)
$\displaystyle=-\bm{\eta}^{T}\frac{\partial\mathbf{M}^{-1}}{\partial
q_{j}}\frac{\partial\operatorname{ID}}{\partial\dot{q}_{i}}-\bm{\eta}^{T}\mathbf{M}^{-1}(\mathbf{q})\frac{\partial^{2}\operatorname{ID}}{\partial\dot{q}_{i}\partial
q_{j}}-\bm{\eta}^{T}\mathbf{M}^{-1}\frac{\partial^{2}\operatorname{ID}}{\partial\dot{q}_{i}\partial{\ddot{\mathbf{q}}}}\frac{\partial\operatorname{FD}}{\partial
q_{j}}$ (8) (Using the noted identity and
$\frac{\partial\operatorname{ID}}{\partial{\ddot{\mathbf{q}}}}=\mathbf{M}$…)
$\displaystyle=\bm{\eta}^{T}\mathbf{M}^{-1}\frac{\partial\mathbf{M}}{\partial
q_{j}}\mathbf{M}^{-1}\frac{\partial\operatorname{ID}}{\partial\dot{q}_{i}}-\bm{\mu}^{T}\frac{\partial^{2}\operatorname{ID}}{\partial\dot{q}_{i}\partial
q_{j}}-\bm{\mu}^{T}\frac{\partial\mathbf{M}}{\partial\dot{q}_{i}}\frac{\partial\operatorname{FD}}{\partial
q_{j}}$ (9) (Using that
$\frac{\partial\operatorname{FD}}{\partial\dot{q}_{i}}=-\mathbf{M}(\mathbf{q})^{-1}\frac{\partial\operatorname{ID}}{\partial\dot{q}_{i}}$
again and that $\frac{\partial\mathbf{M}}{\partial\dot{q}_{i}}=0$…)
$\displaystyle=-\bm{\mu}^{T}\frac{\partial\mathbf{M}}{\partial
q_{j}}\frac{\partial\operatorname{FD}}{\partial\dot{q}_{i}}-\bm{\mu}^{T}\frac{\partial^{2}\operatorname{ID}}{\partial\dot{q}_{i}\partial
q_{j}}$ (10) (Under treating $\bm{\mu}$ as fixed after its initial
evaluation:)
$\displaystyle=-\left(\frac{\partial\operatorname{FD}}{\partial\dot{q}_{i}}\right)^{T}\frac{\partial\mathbf{M}\bm{\mu}}{\partial
q_{j}}-\frac{\partial}{\partial
q_{j}}\left[\bm{\mu}^{T}\frac{\partial\operatorname{ID}}{\partial
q_{i}}\right]$ (11)

We have tried to strike a balance, both here and in the paper, between
presenting the key steps of the derivation, but not overburdening the reader
with it. Regardless, we have double checked the final formulas for accuracy,
and hope that this extra derivation helps bring the reviewer additional
clarity. The CasADi code that implements these calculations is also made
available open source at https://tinyurl.com/468ynkuu for independent
verification.

algorithms that use is approach as”

We thank the reviewer for noting the typo. We have updated the sentence to
read as follows: “algorithms that use this approach as”.

esults: ”implementation in complied C/C++”

We thank the reviewer noting the typo in ‘compiled.’

ould you provide an explanation for Modified RNEA being faster than RNEA only
for $n>=18$ even though its main point was that it needed only 1 RNEA pass
instead of 2.

We note that this was a minor error in implementation whereby we were
incurring unnecessary overhead with the Modified RNEA. (Modified RNEA included
more CasADi function calls, which have high overhead in the MATLAB CasADi
interface. We now ensure that all cases have a single function that interfaces
with CasADi for fairness and accuracy of the comparison.) That the figure has
now been corrected and is attached below as Fig. 10.

Figure 10: The computation time of the partials of the dynamics using the
aforementioned methods.

ould you provide an explanation from where the difference in computation time
between Fig 5 (around 10s to convergence for a swing-up OCP with n = 15) and
Fig 6 (around 1000s for n = 15) comes from? Or be more explicit about the
difference between the two experiments that leads to such a time difference?
Since from what I understand there are both with pendubot.

We appreciate the reviewer bringing this issue to our attention. There was a
discrepancy in the code for these cases. Following the initial review, we have
fixed this discrepancy and re-run all of the results to ensure consistency
across them. Since other reviewers were more interested in performance for low
DoFs, we updated Figure 6 to the case of $n=7$ DoFs for broader appeal.

he last section is interesting as it makes a link with a Kuka system simulated
in Matlab. These implementation efforts toward a real system represent a nice
initiative. However I feel the work done for that is not well exploited in the
paper. The contribution is mainly about reducing the computational cost of DDP
with new variants yet they are all grouped together in Fig 8. We already know
that DDP has quadratic convergence compared to iLQR which has super-linear
convergence (as pointed out in the first sections) so we would get a similar
result if the red curve was plotting ”DDP with tensor contraction” and the
blue curve the classic iLQR method.

Your last part on the simulated Kuka arm should be enhanced to bring
information related to your contributions. You could for instance keep Fig 8
but plot the convergence for each individual variant (to compare them with the
classic methods) or plot instead the cumulated computation time for each
variant to perform the movement of Fig 9 (which would look a bit like Fig 5,
again to compare the new variants and see their advantages for evaluation
time). Or a figure comparing the quality of Fig 9 movement along time with
each variant to show that indeed the new variants keep the benefits of DDP
(second order terms) without sacrifice in evaluation time (as shown in Fig 4
to Fig 7). Due to the nature of your contribution I believe there is no need
to have experimental results on a real system, this study of the simulated
robotic arm is enough if it is properly exploited.

Towards improving the KUKA LBR example, while addressing a main concern raised
by the AE (see AE comment AE.2), we randomize the initial control sequence for
the $7-$link KUKA LBR manipulator to study whether DDP or LQR is more robust
to bad initialization. We present the following findings along with Fig. 11
[2.5em]0em  Since the solution of an OCP is dependent on the initial
conditions, we randomize the initial controller using Ornstein-Uhlenbeck
process noise for the $7-$link KUKA LBR manipulator and solved the OCP problem
using either iLQR/DDP in $40$ separate instances. The first pane of Fig.
LABEL:fig:RandomCtrl illustrates the progression of that manipulator from an
initial configuration to the balanced upright configuration along an optimal
trajectory. Subplot (a) of Fig. LABEL:fig:RandomCtrl illustrates the empirical
probability density function (pdf) of the number of iterations over those
optimizations. Subplot (b) illustrates the pdf of the log of the final cost of
each optimal solution. As illustrated in subplot (a), in $95\%$ of the
simulations, iLQR had a higher number of iterations. On average, iLQR had
three times as many iterations as DDP. From subplot (b), we note that most of
the simulations regardless of DDP/iLQR converged to similar solution; in fact,
DDP converged to a different solution than iLQR in only three instances.
Figure LABEL:fig:RandomCtrl illustrates that the inclusion of second-order
information in DDP will result in similar converged solution and in fewer
iterations. This result is powerful in that the addition of second-order
information results in algorithm that converges in less iterations and to the
same minima.

Figure 11: Simulations on $7-$link KUKA LBR manipulator considering randomized
initial controller given by an Ornstein-Uhlenbeck process noise.

We acknowledge that this comparison does not fully address the role that the
different DDP variants play here, but with the many pieces of feedback to
balance, this approach seemed like the best compromise to us.

We would like to note however, that this additional study showed that iLQR and
DDP almost always converge to the same solution (for the cases considered),
and all the DDP variants take the same steps from iteration to iteration
(since they are based on the same derivatives). In this sense, the main
difference between the variants is their derivative computation time, which is
fully detailed in Figure 4 in the manuscript.

have pointed out the most obvious typos, an additional thorough reading of the
paper could eliminate the other ones if there is any remaining.

We thank the reviewer for noting typos in the paper. We have rectified those
typos and done our best to ensure that the revised paper does not contain any
more.