# Optimizing Complex Control Systems with Differentiable Simulators: A Hybrid Approach to Reinforcement Learning and Trajectory Planning

Amit Parag, Nicolas Mansard, Ekrem Misimi

► To cite this version:

Amit Parag, Nicolas Mansard, Ekrem Misimi. Optimizing Complex Control Systems with Differentiable Simulators: A Hybrid Approach to Reinforcement Learning and Trajectory Planning. 2025 IEEE International Conference on Robotics and Automation (ICRA), IEEE, May 2025, Atlanta, United States. pp.9214-9220, 10.1109/ICRA55743.2025.11127746 . hal-05294931

HAL Id: hal-05294931

https://hal.science/hal-05294931v1

Submitted on 2 Oct 2025

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>
# Optimizing Complex Control Systems with Differentiable Simulators: A
Hybrid Approach to Reinforcement Learning and Trajectory Planning

Amit Parag<fnref n="1" />, Nicolas Mansard<fnref n="2" />, and Ekrem Misimi<fnref n="1" />

**Abstract**— Deep reinforcement learning (RL) often relies on simulators as abstract oracles to model interactions within complex environments. While differentiable simulators have recently emerged for multi-body robotic systems, they remain underutilized, despite their potential to provide richer information. This underutilization, coupled with the high computational cost of exploration-exploitation in high-dimensional state spaces, limits the practical application of RL in the real-world. We propose a method that integrates learning with differentiable simulators to enhance the efficiency of exploration-exploitation. Our approach learns value functions, state trajectories, and control policies from locally optimal runs of a model-based trajectory optimizer. The learned value function acts as a proxy to shorten the preview horizon, while approximated state and control policies guide the trajectory optimization. We benchmark our algorithm on three classical control problems and a torque-controlled 7 degree-of-freedom robot manipulator arm, demonstrating faster convergence and a more efficient symbiotic relationship between learning and simulation for end-to-end training of complex, poly-articulated systems.

## I. INTRODUCTION

In this paper, we present a method that combines reinforcement learning with a differentiable trajectory optimizer to efficiently compute optimal control policies by leveraging both data-driven learning and model-based planning techniques.

The underlying challenge in robotic control is typically framed as an optimal control problem (OCP), modeled as a Markov decision process (MDP), where the objective is to find a policy that maximizes a cumulative reward over a given horizon. There are two main, arguably complementary, approaches that solve this problem: planning-based methods that minimize a cost function over future horizons, and learning-based methods that optimize policies to maximize rewards.

**Planning-based methods**, such as trajectory optimization, compute control sequences by minimizing a cost function over a future horizon. These methods use differentiable simulators (e.g., MuJoCo [1]) to handle non-linearities, dynamics, and constraints [2], and employ techniques such as analytical [3] or automatic differentiation [4], [5]. While effective, these methods can be computationally expensive, especially with algorithms like Differential Dynamic Programming (DDP) [6], [7], which often require many iterations to converge.

**Learning-based methods**, particularly RL, address the control problem by iteratively optimizing a policy to maximize a reward function through interaction with the environ-

ment [8]. Although flexible and adaptable, RL approaches suffer from high sample complexity and sensitivity to hyper-parameters [9]–[11]. Furthermore, RL relies on approximate value function estimates, which are rarely accurate and provide limited guidance for policy optimization [12], [13].

When used in isolation, both approaches have notable drawbacks. Planning-based methods, while precise, are computationally intensive and less effective under model inaccuracies or disturbances. Conversely, RL methods require extensive data and hyperparameter tuning, making them impractical in scenarios where data collection is costly [14], [15].

To address these limitations, several hybrid approaches have been developed to combine the strengths of planning and learning. For example, trajectory optimizers have been used to guide RL agents’ exploration strategies [16]–[18], while learning-based techniques have enhanced trajectory optimization by learning critical quantities like dynamic models or cost functions [19]–[22]. Recent studies have also integrated Sobolev regression [23], [24] to improve learning efficiency and robustness [25]–[27]. These hybrid methods leverage both model-based insights and data-driven adaptability, achieving more sample-efficient and robust solutions.

Our proposed method adopts this hybrid philosophy by coupling a differentiable trajectory optimizer with an RL loop in an actor-critic framework. The RL loop serves as the *actor* learning the policy, while the trajectory optimizer plays the *critic* refining the policy by providing accurate state-value pairs and state-control trajectories. This iterative process is guided by the recursive application of Bellman's optimality principle [28], [29]. By integrating Differential Dynamic Programming (DDP) [6], [30] and employing Sobolev regression [23], [31], our method achieves high accuracy with fewer rollouts, reducing dependency on hyperparameters and improving robustness.

Our contributions are as follows:

* We propose a novel algorithm, Differential Policy Value Programming ($\partial$PVP), which iteratively refines estimates of the global value function and state-control trajectories. $\partial$PVP uses a value function approximator as a proxy for terminal cost in DDP while state and control trajectory approximators are used to warm-starting DDP.

* We demonstrate that $\partial$PVP achieves faster convergence, reduces data requirements and sensitivity to hyperparameters, thereby improving scalability to high-dimensional and complex robotic control tasks.

Using the value function as the terminal cost in DDP

<footnote n="1">¹SINTEF Ocean, Norway, LAAS-CNRS, France</footnote>

<footnote n="2">²LAAS-CNRS, France</footnote>
ensures optimal performance, provided the trajectory ends in a region where the terminal cost accurately represents the value function [21]. This method effectively transforms the problem into an infinite-horizon MPC framework [32], [33].

Bellman's equation generates a time-dependent value function, representing the optimal cost-to-go. This time dependency means that, toward the end of the trajectory, the planning horizon for states becomes very brief, which may lead to short-sighted behavior. To counteract this, a carefully designed terminal cost function, $l^f$, can inform the controller about future events beyond its planning horizon. However, crafting such a custom $l^f$ is challenging. Alternatively, using an approximation of the value function, $V^*$, provides a quasi-infinite prediction horizon, ensuring convergence and offering a significant improvement over traditional methods such as moving horizon techniques [34], receding horizon approaches [35], or those incorporating terminal inequality constraints [36].

The rest of the paper is organized as follows: Section II reviews the foundational concepts of optimal control, Differential Dynamic Programming and Sobolev Regression, Section III describes the $\partial$PVP algorithm, and Section IV presents the evaluations. We conclude in Sec V.

## II. PRELIMINARY

### A. Optimal Control

We consider a time-discrete Markov decision process (MDP) for a system with autonomous dynamics.

In an $n$-dimensional vector space, the transition from the current state $x$ to the next state $x_+$, under control $u$ and environment parameters $\Omega$, is given by:

$$x_+ = f(x, u, \Omega) \quad (1)$$

where $x, x_+ \in \mathcal{X}$ and $u \in \mathcal{U}$.

The objective is to find the state-control pair $(x(t), u(t))$ that minimizes a cost functional $L(X, U)$. For a finite time horizon $T$, this cost functional is expressed as:

$$L(X, U) = \sum_{k=0}^{T-1} l(x_k, u_k) + l^f(x_T) \quad (2)$$

where $l(x, u)$ is the running cost, and $l^f(x_T)$ is the terminal cost. The optimal solution over this horizon is denoted by $(X^*, U^*)$.

### B. Recursive Optimality

The value function $V(x)$ represents the optimal cost starting from state $x$ and following the optimal trajectory $X^*$ with policy $\pi: x \to \pi(x) \in \mathcal{U}$.

The cost-to-go function $V_{T-k}(x)$ over the remaining horizon $T-k$ is given by:

$$V_{T-k}(x) = \min_{U} \sum_{j=k}^{T-1} l(x_j, u_j) + l^f(x_T) \quad (3)$$

Replacing $l^f$ with the value function results in the recursive form:

$$V_{T-k}(x) = \min_{u} [l(x, u) + V_{T-k-1}(f(x, u))] \quad (4)$$

With an accurate estimate of the value function, the problem can be solved as an infinite horizon problem, where the cost-to-go becomes independent of timestep: $V_T(x) = V(x) \quad \forall T > 0$.

### C. Differential Dynamic Programming

Differential Dynamic Programming (DDP) [6], [37] is a second-order iterative algorithm that leverages Bellman's optimality principle to solve nonlinear optimal control problems. It optimizes control sequences by alternately performing a backward and forward pass on the current estimates of state-control trajectories $(X, U)$. In the backward pass, DDP constructs a quadratic approximation of the value function around a nominal trajectory using a Taylor expansion [7], retaining first- and second-order derivatives while discarding higher-order terms. This quadratic form captures the local dynamics of the value function, enabling the algorithm to compute the optimal change in control inputs to minimize the cost-to-go. The forward pass then refines the state-control trajectory by updating the control policy and simulating the system dynamics forward in time.

Warm starts, which initialize the optimization with a well-chosen trajectory, significantly enhance DDP's performance by reducing the number of iterations required for convergence. For instance, in [38], warm starts were used to initialize DDP with trajectories derived from simpler control laws, leading to faster convergence in robotic control tasks. In [39], warm starts were used in conjunction with DDP to handle high-dimensional control systems in real-time control scenarios, improving both stability and convergence speed. These examples show how warm starts are crucial for applying DDP efficiently in complex, high-dimensional control problems where finding an optimal solution from scratch would be computationally prohibitive.

### D. Sobolev Regression

Sobolev regression extends traditional regression by incorporating the derivatives of the function. The goal is to minimize a loss function that includes penalties for deviations in both the function values and their derivatives.

The loss function for Sobolev regression is composed of two parts:

1. 0th Order Loss: This loss measures the discrepancy between the model's predictions and the target values. It is given by:

$$L_0 = \frac{1}{N} \sum_{i=1}^{N} \|y_i - \hat{y}_i\|^2 \quad (5)$$

where $y_i$ represents the true target values, $\hat{y}_i$ represents the predicted values from the model, and $N$ is the number of data points.

2. 1st Order Loss: This loss penalizes the difference between the derivatives of the model's predictions and the target derivatives. It is given by:

$$L_1 = \frac{1}{N} \sum_{i=1}^{N} \|\nabla y_i - \nabla \hat{y}_i\|^2 \quad (6)$$
where $\nabla y_i$ and $\nabla \hat{y}_i$ are the derivatives of the target values and predictions, respectively.

The overall Sobolev loss function combines these two losses:

$$L_{Sobolev} = L_0 + \lambda_d L_1 \quad (7)$$

where $\lambda_d$ is a regularization parameter that balances the contribution of the 1st order loss relative to the 0th order loss.

In our implementation, we focus on the 1st order derivatives with norm $d = 2$ to maintain a manageable computational load. This choice balances accuracy and efficiency by limiting the complexity of the automatic differentiation requirements.

### III. DIFFERENTIAL POLICY VALUE PROGRAMMING

Differential Policy Value Programming (DPVP) integrates Sobolev learning, Bellman's optimality principle, and Differential Dynamic Programming (DDP) to iteratively refine estimates of the global value function and state-control trajectories. The algorithm employs three neural networks:

* A **value network** $V_\alpha^n$ approximates the value function $V(x)$. Given a state $x$, $V_\alpha^n$ predicts the expected cost-to-go from that state, reflecting the cumulative cost from the current state to the end of the planning horizon.

* A **state network** $X_\beta^n$ learns to predict future state trajectories. For a given state $x$, $X_\beta^n$ forecasts the subsequent states $x_{k+1}$ over the planning horizon.

* A **control network** $U_\gamma^n$ estimates the control inputs required at each step. Given a state $x$, $U_\gamma^n$ provides all the subsequent control actions over the planning horizon.

Here, $\alpha, \beta$, and $\gamma$ denote the parameters of the respective networks at the $n$th iteration. Each network takes a single state $x$ as input.

Initially, $\partial$PVP generates a batch of trajectories over a predefined horizon length $T$ without a terminal cost model ($l^f = 0$). These trajectories are used to train $V_\alpha^{n=1}$, $X_\beta^{n=1}$, and $U_\gamma^{n=1}$. Sobolev regression is employed to approximate the value function $V(x)$, while classical regression is used to learn state and control trajectories.

#### A. Iterative Learning

In each iteration, the network $V_\alpha^n$ is updated to act as a proxy for the terminal cost function. The updated cost-to-go function is:

$$V^n(x) = \min_u \left[ \sum_{k=0}^{T-1} l(x_k, u_k) + V_\alpha^{n-1}(x_T) \right] \quad (8)$$

Here, $V_\alpha^{n-1}$ represents the value function from the previous iteration. The state network $X_\beta^n$ and the control network $U_\gamma^n$ provide candidate trajectories during the forward pass in DDP.

As the number of iterations $n$ increases, $V_\alpha^n$ is trained to approximate the value function over a horizon of $T \times n$. After sufficient iterations, $V_\alpha^n$ used as a proxy for the terminal cost will approximate the value function for an infinite horizon

at the terminal step. This refinement stabilizes the solution by ensuring that the terminal cost function more accurately reflects the long-term behavior of the system.

a) **Warm-Starting:** Warmstarting significantly enhances the efficiency of $\partial$PVP. By initializing DDP with trajectories predicted by the neural networks $X_\beta^n$ and $U_\gamma^n$, the algorithm reduces the number of iterations required for DDP to convergence. This approach uses previously learned trajectories to provide a more informed starting point for the optimization, leading to faster and more stable convergence, especially in high-dimensional control problems.

b) **Value function approximation:** The convergence of $V_\alpha^n$ to the infinite horizon value function enhances the stability and robustness of the solution. It reduces the sensitivity of the solution to the choice of the horizon length $T$ and improves the overall performance by providing a more accurate representation of the terminal cost. Consequently, the iterative process not only refines the estimates of the value function but also leads to faster convergence and better scalability for complex control problems.

The pseudocode for $\partial$PVP is shown in Algorithm 1, and outlines the process of sampling trajectories, training neural networks, updating the terminal cost, and warm-starting DDP with the learned trajectories.

#### Algorithm 1: $\partial$PVP

**Algorithm parameters:** Horizon length $T$, iterations $N$, sample size $s$;

**Initialize:** Value network $V_\alpha^0$, state network $X_\beta^0$, control network $U_\gamma^0$;

**Initialize:** DDP solver;

**for** $n = 1$ to $N$ **do**

Sample $s$ trajectories using DDP;

Train $V_\alpha^n$ using Sobolev regression on the sampled data;

Train $X_\beta^n$ using classical regression on the sampled state trajectories;

Train $U_\gamma^n$ using classical regression on the sampled control trajectories;

Update DDP terminal cost: $l^f(x_T) \leftarrow V_\alpha^n(x_T)$;

Warm-start DDP with the updated networks $X_\beta^n$ and $U_\gamma^n$;

#### B. Architecture of the trainable models

$V_\alpha^n$, $X_\beta^n$, and $U_\gamma^n$ were initially designed using a three-tailed feedforward network with shared hidden layers and different output tails. This approach was used to evaluate the trade-off between the longer training times of multi-tailed networks and the potential benefits of shared hidden layers. However, it resulted in increased training time. The practical benefits of shorter training time far outweighed the theoretical advantages of common hidden layers. We do not include those results in this paper. For our experiments, we decided to learn the three different quantities - value function,
state trajectory and control trajectory - on three different feed forward networks.

a) State-Control Approximators: The state and control estimators are implemented with simple deep feed forward networks with 6 hidden layers with ReLU, ELU, Tanh alternately applied. The output and input shape depends on the dimensions of the optimal control problem.

b) Value Function Approximator: Differential Dynamic Programming, as described in Sec II-C, is a second-order algorithm that computes the first-order, $V'_{\alpha}(x)$, and second-order, $V''_{\alpha}(x)$, derivatives of the value function during its backward pass. These derivatives are necessary for refining the state-control trajectories iteratively. However, directly computing $V''_{\alpha}(x)$ for high-dimensional systems, such as a manipulator arm, becomes computationally infeasible.

To address this, we use the Gauss-Newton approximation [40], [41] to simplify the calculation of the value function's second-order derivatives. Specifically, we approximate the value function $V_{\alpha}(x)$ as the squared sum of residuals $R_{\alpha}(x)$:

$$V_{\alpha}(x) = R_{\alpha}(x)^2. \qquad (9)$$

Here, $R_{\alpha}(x)$ is the residual function. By using this formulation, we can express the first and second-order derivatives of the value function as:

$$ \begin{aligned} V'_{\alpha}(x) &= 2R'_{\alpha}(x)^T R_{\alpha}(x), \\ V''_{\alpha}(x) &\approx 2R'_{\alpha}(x)^T R'_{\alpha}(x). \end{aligned} \qquad (10) $$

This approximation simplifies the computation by approximating the Hessian matrix of the value function $V''_{\alpha}(x)$ using only the first-order derivatives of the residual function $R'_{\alpha}(x)$. This method reduces computational complexity while still capturing the essential curvature information necessary for DDP to converge efficiently in high-dimensional spaces.

The update policy of this network is defined in Sec II-D.

## IV. EVALUATION

We benchmark our algorithm on three classical control problems: unicycle [42], cartpole [14], [43] and inverted pendulum [43], along with a torque-controlled 7-dof manipulator arm<fnref n="1" />. Since the precision of our algorithm increases with $(n + 1) \times T$, where n, T are the iterations and predefined horizon of the problem, the quality of the learned value function should depend on the $n^{th}$ iteration or on T or both. We establish validation dataset for simple classical control systems by sampling a large collection of long locally optimal trajectories of horizon length 1000. For the 7-dof manipulator, sampling a similar validation dataset is not feasible.

In this section, we present certain results of interest for unicycle and 7-dof manipulator arm for the complete algorithm<fnref n="2" />. The source code is available at https://github.com/amitparag/DpVp

<figure><img src="image_2.png" /><figcaption>Fig. 1: Comparison of predicted control trajectory and corresponding locally long optimal trajectory, for the manipulator task</figcaption></figure>

<figure><img src="image_3.png" /><figcaption>Fig. 2: Comparison of predicted control trajectory and corresponding locally long optimal trajectory, for the manipulator task</figcaption></figure>

### A. Manipulator Arm

The goal for the manipulator arm is to reach a static target with its End-Effector (EE): the OCP is formulated as a static end-effector pose reaching task, controlled in torque with additional regularization on both state and torque controls. The cost functional also penalizes deviation from state limits. The dynamics of this 7-dof manipulator arm were computed through [3]. We tested the policy trails using the Bullet simulator [44]. During training, we drew 150 locally optimal samples with a horizon length of 200 from the 14-dimensional configuration space in each of the 20 iterations. This resulted in our DDP computing 10347 rollouts overall. The training phase lasted 61 minutes.

Fig. 1 and Fig. 2 show the predicted state trajectories and control trajectories respectively, as compared with a corresponding long, locally optimal trajectory. The advantage of learning state-control in a supervised setting can be seen in the smoothness of the predictions. The predicted state trajectories are nearly coincident with the corresponding long optimal trajectories. When these learned trajectories are used as a reference guide to direct the policy search in DDP, the number of rollouts/attempts required by DDP reduces, as opposed to the number of rollouts required by DDP when not provided any intelligent guesses.

In Fig. 3, we illustrate the number of rollouts required by our converged networks alongside with DDP to solve a problem. We use our trained networks inside our trajec-

<footnote n="1">¹Kuka lwr iiwa R820 14</footnote>

<footnote n="2">²The experiments were performed on core i9 processor with 32 Gb RAM</footnote>
<figure><img src="image_4.png" /><figcaption>Fig. 3: Histogram of attempts (rollouts) required by DDP to perform a task vs attempts required by DDP when guided by learning ($\partial$PVP), shown here for the static pose reaching task with a manipulator arm.</figcaption></figure>

<figure><img src="image_5.png" /><figcaption>Fig. 4: Bellman Residuals for the end effector pose reaching task.</figcaption></figure>

tory optimizer to solve for 600 uniformly sampled initial configurations from the configuration space. For those very same initial configurations, we also use trajectory optimizer but without any help from the trained networks. We then compare the number of rollouts needed in both cases. We see that when learning is used in DDP, the number of attempts made to find optimal policies decrease drastically. The average rollouts required by $\partial$PVP is 2.71, while the average rollouts required by DDP is 9.3.

Fig. 4 shows the difference between the value functions predicted in the $n^{th}$ and $(n+1)^{th}$ iteration. As we see, the relative bellman residuals stabilizes after the $10^{th}$ iteration. This effectively implies that 10 iterations are sufficient for convergence.

## B. Unicycle

Fig. 5 shows the evolution of the predicted value function for the unicycle as $n$ increases. The unicycle, a 3 DOF point mass in the 2D Cartesian horizontal plane, can move forward or rotate with unconstrained longitudinal and angular velocity. The state space is a 3D Cartesian vector space $s = [q_0, q_1, q_2] = [x, y, \theta]^T$ and the control configuration is a two-dimensional vector: $u = [v, \omega]$. The task is to reach the state configuration $q = [0, 0, 0]^T$ while minimizing the weighted cost functional: $L(X, U) = ||w_1q||^2 + ||w_2u||^2$.

<figure><img src="image_6.png" /><figcaption>Fig. 5: Comparison of the topology of value functions predicted by residual network at iterations 1, 5, 10 with value function at horizon length 1000: $V(s, T = 1000)$. The target position is shown in red. The unicycle can start anywhere in the configuration space and tries to reach for target.</figcaption></figure>

To compare the differences in the predictions of the value function, we compute a *quasi-infinite* horizon value function

<figure><img src="image_7.png" /><figcaption>Fig. 6: Errors between the predicted value function at iterations 1, 5, 10 and $V(s, T = 1000)$ for unicycle. The red dot is the target.</figcaption></figure>

<figure><img src="image_8.png" /><figcaption>Fig. 7: Trajectories followed by unicycle under $\partial$PVP and DDP.</figcaption></figure>

validation dataset by sampling for locally long optimal horizons of length 1000. We use this dataset, $V(s, T = 1000)$, to establish distance to infinite horizon in value learning.

As we can see in Fig. 5, the algorithm quickly learns the overall topology of the value function across the state space. As the number of iterations increases, it seems to refine the inherent symmetry in the topology. Fig. 6 shows the difference between the quasi-infinite horizon value function and the predicted value functions at iterations 1, 5, and 9. We observe that the iterative aspect of $\partial$PVP allows it to learn the value function over long horizons quite well despite initialization in a short horizon. However, there are regions in the configuration space that are difficult to handle (shown in white). We suspect that the non-holonomic constraints in the unicycle environment lead to singularities that destabilize learning. With more learning, the algorithm overcomes the presence of singularities. By the $10^{th}$ iteration, the algorithm identifies singularities that are symmetrically distributed around the goal and aligned with the $q_1$ axis in Fig. 6.

## C. Stability of Infinite Horizon

$V_\alpha^n$ approximates value function at different horizons during its iterations. As noted before, the horizon approximated during an iteration is a function of $n$ and $T$. In Fig. 7, we compute the solution to an identical task (for unicycle) over different horizon lengths. We see that the trajectories computed by a converged $V_\alpha^n$ are much closer to the quasi steady-state solution even when the preview horizon is short. Comparatively, without a terminal proxy, the variance of trajectories computed by DDP is much higher.

Fig. 8 shows Bellman residuals and validation errors for systems where it is computationally cheap to generate large validation datasets. The DDP solver computed locally
<figure><img src="image_9.png" /><figcaption>Fig. 8: Bellman residuals and validation errors for systems where it is computationally cheap to generate large validation datasets. The DDP solver computed locally optimal trajectories with $T = 50$ that served as datasets in every iteration.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>Fig. 9: Noise is injected into the system at the initial step. Even under incorrect initialization, Sobolev regression and Bellman's optimality principle force the value estimation to converge.</figcaption></figure>

optimal trajectories of horizon length 50 that served as dataset in every iteration.

Fig. 9 establishes robustness under noise. Noise is injected into the system at the initial step. We observe that even under incorrect initialization, Sobolev regression and Bellman's optimality principle force the value estimation to converge.

In Fig. 10, we see that learning converges faster when $T$ increases. For $T > 40$, convergence is quickly reached. For $T < 40$, convergence is slower.

V. DISCUSSION AND CONCLUSION

In this work, we showed that reformulating reinforcement learning as a combination of iterative supervised learning phase, with emphasis on value functions, and simulators allows for reduction in trials needed to find an optimal solution. The iterative supervised learning phase enforces the stability of predictions through Sobolev regression while learning guided by recursive optimality reduces the dependence on the hyper parameters.

There are a few points of interest in the way the information from policy trials are used in the supervised learning phase. First, in the learning phase the trajectory optimizer plays the role of an oracle that provides an estimation of time dependent value function at each state for every state trajectory, gradients and hessians of each value function along the state trajectory, control trajectory and first order derivatives of control. However, we choose to keep the state with the highest preview horizon for learning the value function. This can be seen as an under-utilization of the state-value information. On the one hand, states with value function estimations of shorter preview horizon can augment

<figure><img src="image_11.png" /><figcaption>Fig. 10: Learning convergence rates as a function of horizon length $T$. For $T > 40$, convergence is reached quickly, while for $T < 40$, convergence is slower.</figcaption></figure>

the dataset but on the other hand, this approach can lead to a noisier dataset which can effectively slow down convergence.

To enforce a shorter training time, we also discarded the use of second order derivatives of value function. This also precluded the use of first-order but high-dimensional derivatives of the control trajectory.

Similarly, we tested the feasibility of adding a kernel loss function based on Bellman's contraction operator on the learning of state and control trajectories. A kernel loss function that constrains trajectories to satisfy the Hamilton-Jacobi-Bellman criteria of optimal sub-structures: sub-solutions of an indefinite horizon optimal control problem should also be optimal solutions to the corresponding definite horizon sub problems, should in theory, elevate the quality of the predicted state and control trajectories. While theoretically beneficial, it showed little improvement while increasing computation time. We leave it to a future work to explore these avenues.

The training time can be further reduced with the observation that the policies learned by the neural networks serve as a guide for DDP. In our experiments, we observe that the precision of learned state and control trajectories increases only slightly as $\partial PVP$ iterates. This allows flexibility in defining the number of training epochs in every iteration, which we need if the dimensionality of the problem under consideration increases. For experiments with manipulator arm, we choose to learn and improve state-control trajectories only in the initial and final iterations. We found this training strategy, where estimations of value functions are refined at every iteration, whereas estimation of state-control trajectories are refined only at start and end to be good enough for our purposes.

In this paper, we presented the foundational aspects of our algorithm. We believe that coupling approach presented in this paper is mature enough to be tested on legged robots such as quadrupeds and mini cheetah in real time. We plan to present those results in a future work.

REFERENCES

[1] E. Todorov, T. Erez, and Y. Tassa, “Mujoco: A physics engine for model-based control,” in *IEEE International Conference on Intelligent Robots and Systems*. IROS, 2012.
[2] M. Huba, S. Skogestad, M. Fikar, M. Hovd, T. Johansen, and B. Rohal'-Ilkiv, "Selected topics on constrained and nonlinear control," *Slovakia, ROSA. Dolný Kubín*, 2011.

[3] J. Carpentier, G. Saurel, G. Buondonno, J. Mirabel, F. Lamiraux, O. Stasse, and N. Mansard, "The pinocchio c++ library: A fast and flexible implementation of rigid body dynamics algorithms and their analytical derivatives," in *IEEE International Symposium on System Integration (SIH)*. SICE, 2019.

[4] Y. Hu, L. Anderson, T.-M. Li, Q. Sun, N. Carr, J. Ragan-Kelley, and F. Durand, "Diffitaichi: Differentiable programming for physical simulation," *ICLR*, 2020.

[5] R. Tedrake et al., "Drake: A planning, control, and analysis toolbox for nonlinear dynamical systems," 2014.

[6] D. Q. Mayne, "Differential dynamic programming—a unified approach to the optimization of dynamic systems," in *Control and Dynamic Systems*. Elsevier, 1973, vol. 10.

[7] Z. Xie, C. K. Liu, and K. Hauser, "Differential dynamic programming with nonlinear constraints," in *2017 IEEE International Conference on Robotics and Automation (ICRA)*. IEEE, 2017, pp. 695–702.

[8] V. Mnih, K. Kavukcuoglu, D. Silver, A. A. Rusu, J. Veness, M. G. Bellemare, A. Graves, M. Riedmiller, A. K. Fidjeland, G. Ostrovski et al., "Human-level control through deep reinforcement learning," *Nature*, vol. 518, no. 7540, 2015.

[9] J. Tan, T. Zhang, E. Coumans, A. Iscen, Y. Bai, D. Hafner, S. Bohez, and V. Vanhoucke, "Sim-to-real: Learning agile locomotion for quadruped robots," in *Proceedings of Robotics: Science and Systems*, Pittsburgh, Pennsylvania, June 2018.

[10] T. Haarnoja, S. Ha, A. Zhou, J. Tan, G. Tucker, and S. Levine, "Learning to walk via deep reinforcement learning," in *Robotics: Science and Systems XV*. RSS, 2019.

[11] T. Haarnoja, A. Zhou, P. Abbeel, and S. Levine, "Soft actor-critic: Off-policy maximum entropy deep reinforcement learning with a stochastic actor," in *ICML*, 2018.

[12] A. Ilyas, L. Engstrom, S. Santurkar, D. Tsipras, F. Janoos, L. Rudolph, and A. Madry, "A closer look at deep policy gradients," in *International Conference on Learning Representations*. ICLR, 2020.

[13] I. Grondman, L. Busoniu, G. A. Lopes, and R. Babuska, "A survey of actor-critic reinforcement learning: Standard and natural policy gradients," *IEEE Transactions on Systems, Man, and Cybernetics, Part C (Applications and Reviews)*, vol. 42, no. 6, pp. 1291–1307, 2012.

[14] R. S. Sutton and A. G. Barto, *Reinforcement learning: An introduction*. MIT press, 2018.

[15] J. Schulman, F. Wolski, P. Dhariwal, A. Radford, and O. Klimov, "Proximal policy optimization algorithms," *arXiv preprint arXiv:1707.06347*, 2017.

[16] S. Levine and V. Koltun, "Guided policy search," in *International Conference on Machine Learning*. ICML, 2013.

[17] N. Mansard, A. DelPrete, M. Geisert, S. Tonneau, and O. Stasse, "Using a memory of motion to efficiently warm-start a nonlinear predictive controller," in *IEEE International Conference on Robotics and Automation*. ICRA, 2018.

[18] K. Lowrey, A. Rajeswaran, S. Kakade, E. Todorov, and I. Mordatch, "Plan online, learn offline: Efficient learning and exploration via model-based control," in *International Conference on Learning Representations*. ICLR, 2019.

[19] I. Lenz, R. A. Knepper, and A. Saxena, "Deepmpc: Learning deep latent features for model predictive control," in *Robotics: Science and Systems*. RSS, 2015.

[20] A. Tamar, G. Thomas, T. Zhang, S. Levine, and P. Abbeel, "Learning from the hindsight plan—episodic mpc improvement," in *IEEE International Conference on Robotics and Automation*. ICRA, 2017.

[21] M. Zhong, M. Johnson, Y. Tassa, T. Erez, and E. Todorov, "Value function approximation and model predictive control," in *Symposium on Adaptive Dynamic Programming and Reinforcement Learning*. IEEE, 2013.

[22] A. Parag, S. Kleff, L. Saci, N. Mansard, and O. Stasse, "Value learning from trajectory optimization and sobolev descent: A step toward reinforcement learning with superlinear convergence properties," in *International Conference on Robotics and Automation*, 2022.

[23] W. M. Czarnecki, S. Osindero, M. Jaderberg, G. Świrszcz, and R. Pascanu, "Sobolev training for neural networks," in *Advances in Neural Information Processing Systems*. Neural IPS, 2017.

[24] K. Hornik, "Approximation capabilities of multilayer feedforward networks," *Neural networks*, vol. 4, no. 2, 1991.

[25] T. M. Mitchell, S. B. Thrun et al., "Explanation-based neural network learning for robot control," in *Advances in Neural Information Processing Systems*. Neural IPS, 1993.

[26] J.-W. Lee and J.-H. Oh, "Hybrid learning of mapping and its jacobian in multilayer neural networks," *Neural computation*, vol. 9, no. 5, 1997.

[27] J. B. Witkoskie and D. J. Doren, "Neural network models of potential energy surfaces: Prototypical examples," *Journal of chemical theory and computation*, vol. 1, no. 1, 2005.

[28] R. Bellman, "Dynamic programming," *Science*, vol. 153, no. 3731, 1966.

[29] S. Peng, "A generalized dynamic programming principle and hamilton-jacobi-bellman equation," *Stochastics: An International Journal of Probability and Stochastic Processes*, vol. 38, no. 2, pp. 119–134, 1992.

[30] C. Mastalli, R. Budhiraja, W. Merkt, G. Saurel, B. Hammoud, M. Naveau, J. Carpentier, L. Righetti, S. Vijayakumar, and N. Mansard, "Crocoddyl: An efficient and versatile framework for multi-contact optimal control," in *International Conference on Robotics and Automation*. ICRA, 2020.

[31] H. Son, J. W. Jang, W. J. Han, and H. J. Hwang, "Sobolev training for the neural network solutions of pdes," *arXiv preprint arXiv:2101.08932*, 2021.

[32] H. Chen and F. Allgöwer, "A quasi-infinite horizon nonlinear model predictive control scheme with guaranteed stability," *Automatica*, vol. 34, no. 10, pp. 1205–1217, 1998.

[33] B. Hu and A. Linnemann, "Toward infinite-horizon optimality in nonlinear model predictive control," *IEEE Transactions on Automatic Control*, vol. 47, no. 4, pp. 679–682, 2002.

[34] E. Polak and T. Yang, "Moving horizon control of linear systems with input saturation and plant uncertainty part 1. robustness," *International Journal of Control*, vol. 58, no. 3, pp. 613–638, 1993.

[35] D. Mayne and H. Michalska, "An implementable receding horizon controller for stabilization of nonlinear systems," in *29th IEEE Conference on Decision and Control*. IEEE, 1990, pp. 3396–3397.

[36] H. Michalska and D. Q. Mayne, "Robust receding horizon control of constrained nonlinear systems," *IEEE transactions on automatic control*, vol. 38, no. 11, pp. 1623–1633, 1993.

[37] L.-z. Liao and C. A. Shoemaker, "Advantages of differential dynamic programming over newton's method for discrete-time optimal control problems," Cornell University, Tech. Rep., 1992.

[38] Y. Tassa, T. Erez, and E. Todorov, "Synthesis and stabilization of complex behaviors through online trajectory optimization," *2012 IEEE/RSJ International Conference on Intelligent Robots and Systems*, pp. 4906–4913, 2012.

[39] I. Mordatch, E. Todorov, and Z. Popović, "Interactive control of diverse complex characters with neural networks," *Advances in Neural Information Processing Systems*, vol. 28, 2015.

[40] J. Nocedal and S. J. Wright, *Numerical Optimization*, 2nd ed., ser. Springer Series in Operations Research and Financial Engineering. New York, NY, USA: Springer Science & Business Media, 2006, chapter 10 discusses the Gauss-Newton method in detail.

[41] Å. Björck, "Numerical methods for least squares problems," *SIAM Review*, vol. 38, no. 1, pp. 123–125, 1996, comprehensive discussion of numerical methods for solving least-squares problems, including the Gauss-Newton method.

[42] S. Fleury, P. Soueres, J.-P. Laumond, and R. Chatila, "Primitives for smoothing mobile robot trajectories," *Transactions on Robotics and Automation*, vol. 11, no. 3, 1995.

[43] G. Brockman, V. Cheung, L. Pettersson, J. Schneider, J. Schulman, J. Tang, and W. Zaremba, "Openai gym," 2016.

[44] E. Coumans and Y. Bai, "Pybullet, a python module for physics simulation for games, robotics and machine learning," http://pybullet.org, 2016–2021.