Offline Supervised Learning V.S. Online Direct Policy Optimization: A Comparative Study and A Unified Training Paradigm for Neural Network-Based Optimal Feedback Control

As mentioned in Section 2.3, the objective (3) of direct policy optimization with fully-known dynamics can be optimized through stochastic gradient descent. There are two ways to compute the required (stochastic) gradient: back-propagation and the adjoint method.

In back-propagation [24], we first discretize the continuous trajectory with a given initial state $\bm{x}_{0}$ by an ODE integration scheme and get a discrete trajectory to approximate the total cost $J(\bm{u};\bm{x}_{0})$. Based on discrete trajectories, we apply back-propagation through the operations of the trajectories to compute the gradient with respect to the parameters $\bm{\theta}$. Adaptive ODE solvers, which automatically find suitable stepsize under certain error tolerance, can often achieve a better approximation efficiency of the continuous trajectory than a fixed-stepsize integrator. However, directly leveraging adaptive ODE solvers leads to an implicit and unknown depth of the computation graph, resulting in a deterioration of performance and computation time in back-propagation. To overcome this limitation, [38] suggest deleting the computation graph of adaptive solvers and storing the final computed stepsize only, which reduces the memory cost and benefits in the final performance.

The adjoint method [25] takes another avenue to calculate the gradient of the total loss $J$ with respect to parameters $\bm{\theta}$ in a continuous viewpoint. [6] defines the adjoint $\bm{\alpha}(t)=\partial J/\partial\bm{x}(t)$ with an augmented ODE $\frac{d\bm{\alpha}(t)}{dt}=-\bm{\alpha}(t)^{\top}\frac{\partial\bm{f}(\bm{x},\bm{u})}{\partial\bm{x}}-\frac{\partial L(\bm{x},\bm{u})}{\partial\bm{x}}$, where the notations are the same as OCP (1). Then the gradient can be calculated by $\frac{\partial J}{\partial\bm{\theta}}=\int_{0}^{T}\bm{\alpha}(t)^{\top}\frac{\partial\bm{f}}{\partial\bm{u}}\frac{\partial\bm{u}}{\partial\bm{\theta}}+\frac{\partial L}{\partial\bm{u}}\frac{\partial\bm{u}}{\partial\bm{\theta}}dt$, which requires a forward ODE to solve the state $\bm{x}(t),t\in[0,T]$ and a backward ODE with $\bm{\alpha}(T)=\partial M/\partial\bm{x}(T)$ to solve the adjoint $\bm{\alpha}(t)$. The numerical error and memory usage of the adjoint method can be less than the back-propagation; however, the computing time is longer due to the additional calls on ODE solvers. The idea of the adjoint method has been followed by Neural ODE [39] to view the discrete neural network model as a continuous flow. Extending the implementation of [39] back to control problems, efforts are made in neural network-based control under deterministic dynamics [4, 6]. Unfortunately, naively leveraging the adjoint method may cause catastrophic divergence, even in the simple linear quadratic regulator (LQR) case [6], which can be alleviated by adding checkpoints to ODE trajectories [6, 38, 40]. The aforementioned idea of adaptive solvers can also be applied and the stored states can serve as the checkpoints [38].

We compare direct policy optimization with back-propagation or the adjoint method in the optimal landing problem on $\tilde{\mathcal{S}}_{\text{quad}}$ with $T=4,8,16$ in Table 5. The adjoint method costs time about twice as back-propagation. We round up to 4 decimal places to show the slight differences in $T=4,8$. When the time horizon is $T=16$, the adjoint method deteriorates even with checkpoints added. Since memory is not the bottleneck, we use back-propagation with competitive performances and faster training speed to train the direct optimization in the main content.

In this subsection, we compare the performance of direct policy optimization with and without fully-known dynamics. We focus on the optimal attitude control problem of a satellite. With fully-known dynamics, we optimize by back-propagating through the computational graph of the entire trajectory. With unknown dynamics, we utilize the on-policy, model-free method, Proximal Policy Optimization (PPO) [26], as a representative method in reinforcement learning, which has been empirically demonstrated to achieve state-of-the-art results in a wide range of tasks, including robotics [41] and large language models [42].

To ensure a fair comparison, we conduct experiments with RL using either the same number of trajectories as direct policy optimization with fully-known dynamics or five times as many, denoted as 1x and 5x samples, respectively. As presented in Table 6, direct policy optimization with fully-known dynamics has a cost ratio of $1.048\pm 0.034$, which is comparable to optimal control. In contrast, RL (5x Samples) has a worse cost ratio of $1.423\pm 0.434$, and RL (1x Samples) is significantly worse of $12.58\pm 8.959$. These results illustrate that RL algorithms are significantly less sample-efficient.

We demonstrate the trajectories controlled by networks trained using the above methods in Figure 7, starting from the same initial point. While the RL solutions seem plausible, they are less optimal than the ones learned through dynamics. In this case, the optimal cost is 4.03, and cost of direct optimization with fully-known dynamics is 4.06, which are very close. However, the cost of RL (5x) is 5.75 and the cost of RL (1x) is much worse, 27.20. The advantages of using fully-known dynamics are straightforward, as we possess more information about the trajectories. Further details of our experiments are provided in B.1.

In this section, we present robustness experiments in quadrotor’s optimal landing problem, when there is state measurement noise in the close-loop simulation. We only fine-tune 100 iterations within 10 minutes based on the network pre-trained by supervised learning.

Figure 8 and 9 demonstrate that supervised learning outperforms direct policy optimization under moderate disturbances. Additionally, we highlight the superior robustness of the Pre-train and Fine-tune strategy compared to both direct policy optimization and supervised learning. When the noise is small, supervised learning is still able to control the quadrotor and maintain its priority over direct optimization. However, since supervised learning is not designed for stochastic optimal control, it does not suit the case when the noise is very large. Also, the long time horizon enlarges the influence brought by disturbance, accumulates the error along the trajectory, and results in the distribution mismatch phenomenon and poor performances. Fine-tuning brings robustness to the pre-trained network, which can endure much larger disturbance and performs the best even when the performance of supervised learning decreases dramatically as the bottom row in Figure 8 shows. Figure 9 also shows that fine-tuning is able to improve the robustness regardless of the dataset distribution where it is pre-trained on.

In this subsection, we aim to quantify the difficulty of different problems by evaluating the performance of various controllers. We begin by computing the total cost of uncontrolled systems on the test dataset. To provide a benchmark for comparison, we also consider the Linear Quadratic Regulator (LQR) as a baseline. The LQR controller is designed to optimize the system’s performance by linearizing the full dynamics around the terminal state [43]. In the context of satellite altitude control, the terminal point is defined as $\bm{x_{f}}=0$ and $\bm{u_{f}}=0$, while in the quadrotor’s landing problem, it is $\bm{x_{f}}=0$ and $\bm{u_{f}}=mg$. We apply Drake [44] to solve the related finite horizon LQR problems. Following the same problem setup as discussed in the main content, the feedback control $\bm{u}(t,\bm{x})=-K(t)\bm{x}$ is utilized, where $K(t)$ represents a time-varying control matrix, which comes from solving the related continuous-time Riccati Equation. The results of the zero controller, LQR controller, and other controllers trained by different methods (as presented in Table 7) are collectively analyzed for a comprehensive comparison. The performance of both zero control and LQR is inferior, indicating the increased difficulty of the respective problems. Notably, the satellite control problem emerges as the least challenging among the scenarios explored in this study, while the quadrotor problem, especially with a longer time horizon, presents a greater level of difficulty.

Hyper-parameters are summarized in Table 8, and we briefly explain some details.

In our RL experiments, we use a discretization timestep of $\Delta t=0.005$ to build the environment, resulting in 4000 steps per episode. The reward at each timestep $t$ is the running cost integrated over a small time interval, $R_{t}=-\int_{t}^{t+\Delta t}{L}(\bm{x},\bm{u})dt$, and the terminal reward additionally subtracts the terminal cost $M(\bm{x}(T))$. In RL, the goal is to learn the policy $\bm{u}$ that maximizes the expected discounted return $G_{t}=\sum_{k=0}^{\infty}\gamma^{k}R_{t+k+1}$, where the discount factor $\gamma=1.0$ in our finite horizon problem. Note that the expected discounted return in RL to maximize is the total cost in OCP (1) to minimize. We use the PPO algorithm implemented in stable-baselines3 [45].⁵⁵5https://github.com/DLR-RM/stable-baselines3 This algorithm employs the generalized advantage estimation (GAE, 46) technique to reduce variance in the gradient estimate. We experiment with gae_lambda from {1.0, 0.99, 0.98, 0.97, 0.96, 0.95}, ultimately settling on gae_lambda=0.98 as it yields the best results. To ensure fair comparisons, we set the number of inner loop epochs in PPO to 1, so that the collected data is only used once in RL, just as in direct optimization with dynamics. It is worth noting that the number of trajectories used in direct optimization with dynamics is the product of its batch size and number of iterations (1024 $\times$ 2000), which is the same as the number of episodes used in RL (1x) (2048000), as shown in Table 8.

We use Dopri5 (Dormand–Prince Runge–Kutta method of order (4)5) in the satellite problem and RK23 in the quadrotor problem as the ODE solvers. The implementation of direct policy optimization in this paper includes the aforementioned advanced techniques of adaptive solvers and adding checkpoints to achieve desired performance. The adaptive solvers in training are implemented by [38].⁶⁶6https://github.com/juntang-zhuang/torch-ACA Since our main results are architecture-agnostic, we use a fully connected neural network as the closed-loop controller, which takes both time and state as the input, i.e., $\bm{u}^{\mathrm{NN}}(t,\bm{x})$. We use time-marching in satellite and space-marching in quadrotor to successfully generate the BVP solution. We randomly generate a validation dataset of 100 initial states for each problem, and adopt BVP solvers to compute their optimal costs for comparison. As for IVP-enhanced sampling, we follow the same settings as [29] does, with temporal grid points on $0<10<14<16$. The number of initial states in uniform sampling and adaptive sampling is the same. We remark that a learning rate decay strategy is incorporated to achieve good performance in SL, following common practices in the field. However, in the case of direct policy optimization, we are unable to employ similar learning rate decay due to its inefficient optimization process. Given the limited number of iterations currently used, reducing the learning rate in direct policy optimization would significantly slow down the training process and hurt performance. On the other hand, at the fine-tuning stage, we have the flexibility to utilize much smaller learning rates, allowing us to achieve optimal performance.

We present further detailed information on the analysis of optimization landscapes in Figure 3, following the settings outlined in [37]. Our focus is primarily on local evaluations conducted during the training process. We denote the loss function as $l(\theta)$ and the learning rate (used in the supervised learning or direct policy optimization) as $lr$. During the local evaluation, we perform updates similar to stochastic gradient descent (SGD), $\theta^{\prime}=\hat{\theta}-{step\_size}\times\nabla_{\theta}l(\hat{\theta})$, where $step\_size=step\_ratio\times lr$ and $\hat{\theta}$ denotes the current parameter in training. Specifically, the multiple values of $step\_ratio$ selected for supervised learning (SL) are within [1/100, 100], while for direct optimization, the values of $step\_ratio$ are within [1/100, 1/10]. Notably, in direct optimization which is randomly initialized, we are limited to choosing $step\_ratio$ values much smaller than 1.0. This restriction arises due to the differing update rules between Adam (the training optimizer) and SGD (the update rule in local evaluation). The need for this careful adjustment is driven by the high sensitivity of direct policy optimization to the choice of optimizer and learning rate, particularly in the initial stages when gradients are significantly larger. The $step\_ratio$ is not strictly limited to such small values when the weights are pre-trained and have a much smoother landscape. We remark that we use the same $step\_ratio$ in direct optimization for both randomly initialized and with pre-trained weights. By selecting appropriate $step\_ratio$ values within this constrained range, we ensure a valid evaluation process.

The column Data denotes the time for generating the whole dataset by BVP solvers according to the related marching methods. We apply parallel techniques to fasten the total procedure of data generation in quadrotor with 24 processors. The addition time in Quad_Large_Ada is related to the network training in adaptive sampling. The other columns represent the training time of supervised learning, the fine-tuning time based on the pre-trained network, and the training time by direct optimization, respectively. The rows Quad_Small_T4, Quad_Small_T8 and Quad_Small_T16 denote the experiments initiated in the small domain $\tilde{\mathcal{S}}_{\text{quad}}$ with varying time spans of $T=4,8,16$. The last two rows are experiments in $\mathcal{S}_{\text{quad}}$ with $T=16$, where Quad_Large_Uni denotes that the dataset for supervised learning is generated based on uniform initial states, whilst Quad_Large_Ada denotes that the dataset is generated by adaptive sampling.

Table 9 demonstrates that supervised learning costs less time than direct optimization and highlights the priority of the proposed Pre-train and Fine-tune strategy in terms of time consumption. Direct optimization is harder to train and costs more time, especially when the problem is more complex and the time horizon is longer. There is no need for fine-tuning in satellite, since the pre-trained network is already very close to the optimal control. The fine-tuning time reported is only for 100 iterations, which are all within 10 minutes. As shown in Table 3 and Table 4, with only 100 iterations of fine-tuning, the performance of the controller can be improved significantly. We also record the fine-tuning time of 1000 iterations, which is about 92 minutes for both the uniform dataset and the adaptive dataset. All experiments related to fine-tuning are fine-tuned for 100 iterations, except the ones mentioned in Table 4.

We utilize 256 parallel environments to collect data for reinforcement learning. The training process takes 85 minutes in RL (1x) and 434 minutes in RL (5x). The training times are comparable between RL (1x) and BP with fully-known dynamics when using the same number of trajectories, and are approximately five times longer for RL (5x) at a larger sample scale.

We introduce the full dynamics of the satellite problem [3, 33, 11] that are considered in Section 4.2.

The state is $\bm{x}=\left(\bm{v}^{\mathrm{T}},\bm{\omega}^{\mathrm{T}}\right)^{\mathrm{T}}$, where $\bm{v}=(\phi,\theta,\psi)^{\mathrm{T}}$ denotes the attitude of the satellite in Euler angles and $\bm{\omega}=\left(\omega_{1},\omega_{2},\omega_{3}\right)^{\mathrm{T}}$ represents the angular velocity in the body frame. $\phi$ (roll), $\theta$ (pitch), and $\psi$ (yaw) are the rotation angles around the body frame.

The closed-loop control is $\bm{u}(t,\bm{x}):\left[0,T\right]\times\mathbb{R}^{6}\rightarrow\mathbb{R}^{m}$. The output dimension $m$ represents the number of momentum wheels, which is set to 3 as the fully actuated case. Then we define the constant matrix $\bm{B}\in\mathbb{R}^{3\times m}$, $\bm{J}\in\mathbb{R}^{3\times 3}$ to denote the combination of the inertia matrices of the momentum wheels, and $\bm{h}\in\mathbb{R}^{3}$ to denote the rigid body without wheels as

Finally, the full dynamics are

where the skew-symmetric matrix and rotation matrices $\bm{S}(\bm{\omega}),\bm{E}(\bm{v}),\bm{R}(\bm{v}):\mathbb{R}^{3}\rightarrow\mathbb{R}^{3\times 3}$ are defined as

The running cost and the final cost are defined as

where the weights are set as $W_{1}=\frac{1}{2},W_{2}=5,W_{3}=\frac{1}{4},W_{4}=\frac{1}{2},W_{5}=\frac{1}{2}$.

We describe the full dynamics of the optimal landing problem of a quadrotor [27, 29, 34, 35, 36] which are considered in Section 4.3.

The state is $\bm{x}=\left(\bm{p}^{\mathrm{T}},\bm{v}_{b}^{\mathrm{T}},\bm{\eta}^{\mathrm{T}},\bm{w}_{b}^{\mathrm{T}}\right)^{\mathrm{T}}\in\mathbb{R}^{12}$ where $\bm{p}=(x,y,z)^{\mathrm{T}}\in\mathbb{R}^{3}$ denotes the position in Earth-fixed coordinates, $\bm{v}_{b}=(v_{x},v_{y},v_{z})^{\mathrm{T}}\in\mathbb{R}^{3}$ denotes the velocity with respect to the body frame, $\bm{\eta}=(\phi,\theta,\psi)^{\mathrm{T}}\in\mathbb{R}^{3}$ denotes the attitude in Earth-fixed coordinates, and $\bm{w}_{b}\in\mathbb{R}^{3}$ denotes the angular velocity in the body frame.

The closed-loop control is $\bm{u}(t,\bm{x}):\left[0,T\right]\times\mathbb{R}^{12}\rightarrow\mathbb{R}^{4}$. In practice, the individual rotor thrusts $\bm{F}=\left(F_{1},F_{2},F_{3},F_{4}\right)^{\mathrm{T}}$ is applied to steer the quadrotor, where $\bm{u}=E\bm{F}$. We use $l$ as the distance from the rotor to the quadrotor’s center of gravity, and $c$ as a constant that relates the rotor angular momentum to the rotor thrust (normal force), then the matrix $E$ is defined as

With $\bm{F}^{*}=E^{-1}\bm{u}^{*}$, the optimal $\bm{F}^{*}$ can be directly computed, once we obtain the optimal control $\bm{u}^{*}$.

Using the same parameters as [35], the mass and the inertia matrix are defined as $m=2kg$ and $\bm{J}=\operatorname{diag}\left(J_{x},J_{y},J_{z}\right)$, where $J_{x}=J_{y}=\frac{1}{2}J_{z}=1.2416\mathrm{~{}kg}\cdot\mathrm{m}^{2}$. Another constant vector is $\bm{g}=(0,0,g)^{\mathrm{T}}$, where $g=9.81\mathrm{~{}m}/\mathrm{s}^{2}$ denotes the acceleration of gravity on Earth. The rotation matrix $\bm{R}(\bm{\eta})\in SO(3)$ denotes the transformation from the Earth-fixed coordinates to the body-fixed coordinates, and the attitude kinematic matrix $\bm{K}(\bm{\eta})$ relates the time derivate of the attitude with the associated angular rate as

Finally, the dynamics are

where the constant matrices $A$ and $B$ are defined as

The running cost and the final cost are defined as

where $Q_{u}=\operatorname{diag}(1,1,1,1)$, $Q_{pf}=5I_{3}$, $Q_{vf}=10I_{3}$, $Q_{\eta f}=25I_{3}$ and $Q_{wf}=50I_{3}$, where $I_{3}$ denotes the identity matrix in three-dimensional space.