Real-time Neural MPC: Deep Learning Model Predictive Control for Quadrotors and Agile Robotic Platforms

In plain words

A controller that plans ahead, such as model predictive control (MPC), works best with an accurate picture of how the robot moves. Simple physics formulas miss effects like air flowing around a fast drone, but a neural network can learn those effects. The catch is speed: the controller has to re-plan many times per second, and a large network makes every re-plan too slow.

This paper's idea is to keep the network out of the heavy part of the calculation. The controller only needs a straight-line guess of how the network's output changes near the current plan. That guess is cheap to compute for all points of the plan at once, on a CPU or GPU, and the controller then does one quick update per tick instead of solving the whole problem again. The result is that much bigger networks can run inside the control loop in real time. The authors tested this on a fast quadrotor, both in simulation and in flight, and also for a model of the air close to the ground.

TL;DR

• Neural dynamics models are added as a learned residual on top of a first-principles quadrotor model eq. 6 eq. 7.
• The network is linearized at each point of the planned trajectory, and those derivatives are computed in one batch §IV-A §IV-B.
• The controller runs one real-time SQP iteration per tick: prepare the quadratic program while waiting for the sensor, then solve it when the new state arrives §IV-B fig. 2.
• Author-reported: larger models than prior neural MPC at real-time rates, and up to 82% lower positional tracking error than nominal MPC in flight tab. IV.
• Code is an MIT-licensed framework that the authors have since merged into a successor project; the experiment code is GPL-3.0 (see Paper vs code).

Problem

Gradient-based MPC needs a dynamics model inside its optimisation. A simple first-principles model is cheap but leaves out effects such as aerodynamics and friction. A neural model is more accurate but expensive, because the optimiser needs its values and derivatives at every shooting node in every iteration §I §IV-A. Earlier neural MPC work therefore used small networks, Gaussian processes with few support points, or closed the loop without online optimisation §I tab. I.

The paper assumes the learned model only has to be accurate in a neighbourhood of the current iterate, not over the whole training distribution §IV-A. It also assumes that solving the optimisation to convergence is too slow on agile platforms, so a single SQP iteration per control step is used §IV-B.

Method

The method has three parts: a model split, a real-time-iteration (RTI) controller and a batched approximation of the learned part.

Model split. The dynamics are split into a nominal part and a learned part, f = f_F + f_D §III eq. 1. For the quadrotor, the nominal part is a 13-dimensional rigid-body model whose inputs are the four rotor thrusts §VI eq. 6. The learned part adds residual linear and angular accelerations produced by two networks §VI eq. 7. Two simplified residuals are also evaluated: one that depends only on body velocity, and one that adds the commanded inputs §VI eq. 8. A ground-effect residual takes a local patch of a height map, cut around the quadrotor, as additional input §VI.

Optimisation problem. The continuous problem is discretised with direct multiple shooting over N steps of length δt, with RK4 integration (four evaluations of f per step) §III eq. 1. The resulting nonlinear program is solved with sequential quadratic programming (SQP) §III.

Real-time iteration. Each SQP step builds a quadratic program whose continuity rows use the linearised dynamics, and whose consistency rows bound the constraints §IV-A eq. 2 eq. 3. The paper splits one control tick into three phases §IV-B fig. 2:

1. Data-driven preparation: the learned model's value and Jacobian at every node of the current iterate are computed in a batch on CPU or GPU §IV-B.
2. QP preparation: the continuity condition is built from these first-order terms; the Hessian may be approximated by Gauss-Newton under the stated conditions §IV-A §IV-B.
3. Feedback response: when the new measured state arrives, the prepared QP is solved with it as initial state, the iterate is updated, and the first command is sent to the actuators §IV-B.

Approximation. The learned map is replaced near the iterate by a Taylor expansion: a value term, a Jacobian term and, optionally, a Hessian term §IV-A eq. 4. The first-order version is the one used in the RTI scheme §IV-B. The paper states that the induced error is of second order in the change of state and control between nodes for the first-order case, and of third order for the second-order case §IV-A. The authors observe that this error is negligible for agile platforms with small δt, but this is shown by experiment only §IV-A §VII.

Implementation. The authors use CasADi and acados as the optimisation stack and PyTorch for the networks §IV-C. The naive baseline builds the full network inside CasADi as weight matrices and activations and differentiates it there §IV-C.

Equations (as numbered in the compiled paper, see Pitfalls):

\min_{u}\ \sum_{k=0}^{N-1}\mathcal{L}(x_k,u_k)\ \text{ s.t. }\ x_0=\bar{x}_0,\ x_{k+1}=\phi(x_k,u_k,f,\delta t),\ f=f_F+f_D,\ g(x_k,u_k)\le 0 \quad [\text{eq. 1}]

\Delta x_{k+1}=A_k\Delta x_k+B_k\Delta u_k+\bar{\phi}_k-x_{k+1},\quad k=0,\ldots,N-1 \quad [\text{eq. 2}]

-\bar{g}_k\ \ge\ G^x_k\Delta x_k+G^u_k\Delta u_k,\quad k=0,\ldots,N \quad [\text{eq. 3}]

f_{\mathcal{D}}^*(x,u)\approx \bar{f}_{\mathcal{D}}^{\,i}+J^i_{\mathcal{D},k}\begin{bmatrix}x-x^i_k\\u-u^i_k\end{bmatrix}+\tfrac12\begin{bmatrix}x-x^i_k\\u-u^i_k\end{bmatrix}^{\!\top}H^i_{\mathcal{D},k}\begin{bmatrix}x-x^i_k\\u-u^i_k\end{bmatrix} \quad [\text{eq. 4}]

The double-integrator experiment uses f_F = [ṗ; u] with the network output fixed to zero eq. 5 §V.

\dot{x}=\begin{bmatrix}\dot p_{WB}\\\dot q_{WB}\\\dot v_{WB}\\\dot\omega_B\end{bmatrix}=\begin{bmatrix}v_W\\q_{WB}\cdot[0,\ \omega_B/2]^\top\\\tfrac1m\,q_{WB}\odot T_B+g_W\\J^{-1}\left(\tau_B-\omega_B\times J\omega_B\right)\end{bmatrix} \quad [\text{eq. 6}]

f_{\mathcal{D}}=\begin{bmatrix}0_2\\f_{\mathcal{D}_\theta}(x,u)\\f_{\mathcal{D}_\psi}(x,u)\end{bmatrix} \quad [\text{eq. 7}]

f_{\mathcal{D}_a}=\begin{bmatrix}0_2\\f_{\mathcal{D}_\theta}(v_B)\\0\end{bmatrix},\qquad f_{\mathcal{D}_{a,u}}=\begin{bmatrix}0_2\\f_{\mathcal{D}_\theta}(v_B,u)\\0\end{bmatrix} \quad [\text{eq. 8}]

Cost and constraints: the stage cost is quadratic in the state and input deviations from a reference, with inputs bounded between zero and a maximum §VI. Weight matrices and bounds are not given numerically in the paper.

# One control cycle of Real-time Neural MPC  [§IV-B] [fig. 2]
# Input:  iterate omega (state and input per node), network D, sensed state
# Output: first input command u0, updated iterate omega

def control_cycle(omega, dynamics_net, sensed_state):
    # Phase 1, data-driven preparation (batched, CPU or GPU)  [§IV-B]
    # Value and Jacobian of the network at every node  [§IV-A]
    f_bar, jacobian = batched_jacobian(dynamics_net, omega)
    # Phase 2, QP preparation: linearise the continuity rows  [§IV-A] [eq. 2]
    # Uses the first-order terms from phase 1 and RK4 integration  [§III]
    qp = build_qp(omega, f_bar, jacobian)
    # Phase 3, feedback response: solve with the sensed state as x_0  [§IV-B]
    # One SQP iteration only, so the step is a single QP solve  [§IV-B]
    step = solve_qp(qp, initial_state=sensed_state)
    # Update the iterate and return the first input  [§IV-B]
    omega = omega + step
    return first_input(omega), omega

Hyperparameters

Values are as stated in the paper. Where the paper is silent, the row says so.

Results

Author-reported. All numbers below are the paper's, not re-run here.

Setup: simplified quadrotor simulation (single CPU core, network sizes N-L-H, Circle and Lemniscate tracks) §VII-A tab. III; runtime on Jetson Xavier NX, an i7 CPU and an RTX 3000 GPU, two-layer and deeper MLPs and a ResNet-style CNN §V tab. II; real-world flights on the agilicious platform, Circle and Lemniscate tracks, five runs per experiment §VII-B tab. IV.

Runtime, control frequency in Hz (higher is better) tab. II. The "None" column is the no-model baseline per platform. Model columns follow the paper's header: layers and neurons per row, parameter counts, and a CNN.

Simplified simulation, mean Euclidean tracking error in mm tab. III. Column "t" is the average optimisation time in ms. Bold in the source marks the best per column and is not reproduced.

Real-world flights, mean tracking error in mm tab. IV. The paper reports ± standard deviation for some rows. Naive integration crashed in both tracks.

Headline claims, as stated: the abstract reports models over 4000 times larger in parametric capacity than prior neural MPC in a 50 Hz window on an embedded platform, and up to 82% lower positional tracking error than nominal MPC [abstract]. The text gives a factor of about 60 on the embedded CPU and about 4000 on the embedded GPU §V fig. 3. The ground-effect residual reduces vertical tracking error by 72% over a table in one evaluation with 8 flyovers §VII-B fig. 6.

Paper vs code

The official framework repository is TUM-AAS/ml-casadi (MIT, cloned at code_commit). The README marks it as deprecated and says its functionality has been merged into the successor L4CasADi README.md README.md. It is the implementation the paper describes in §IV-C.

The experiment code, TUM-AAS/neural-mpc, is also linked from the paper [abstract]. Its LICENSE is GPL-3.0. Nothing from it is quoted here. It was read via WebFetch, not cloned, so line numbers are not given; cite by file path only.

Findings:

• Early stopping. Paper: early stopping on a validation set §VII. The experiment training script saves the model after every epoch and does not select a best checkpoint or stop early. Its default epoch count is 1000 (neural-mpc ros_dd_mpc/src/model_fitting/mlp_fitting.py, read via WebFetch). Paper and script differ. The paper's stopping rule is the one to follow when reproducing the paper.
• Optimiser and batch size. Paper: Adam, learning rate 1e-4, batch 64 §VII. Script: Adam, learning rate 1e-4 fixed, batch size default 64, no scheduler (neural-mpc ros_dd_mpc/src/model_fitting/mlp_fitting.py, WebFetch). Match.
• Network default. Script default is one hidden layer of 32 units with tanh activation, which is not any of the paper's N-L-H configurations (neural-mpc mlp_fitting.py, WebFetch). Use the paper's table for the experiments.
• Input and output normalisation. Not in the paper. The script wraps the network with mean and standard-deviation normalisation fitted on the training set (neural-mpc mlp_fitting.py, WebFetch).
• Solver. Paper: one SQP step per tick, real-time iteration §IV-B. Code: SQP_RTI is the default; Gauss-Newton Hessian; HPIPM full-condensing QP solver; explicit Runge-Kutta integrator; stage count not confirmed (neural-mpc ros_dd_mpc/src/quad_mpc/quad_3d_optimizer.py, WebFetch). Consistent with the paper, apart from the unconfirmed integrator order.
• Cost. Paper: quadratic stage cost, weights not stated §VI. Code: least-squares cost with diagonal Q (position weights 10, others 0.05 to 0.1) and diagonal R (0.1); terminal cost off by default (neural-mpc quad_3d_optimizer.py, WebFetch). These are code defaults, not paper values.
• Horizon. Paper: N = 10 for the runtime study §V. Code default: N = 20 nodes over 1 s, so 0.05 s steps (neural-mpc quad_3d_optimizer.py, WebFetch). Code default differs from the paper's runtime study.
• Linearisation point. Paper linearises around the current iterate §IV-A §IV-B. Code computes the approximation inputs from the warm-start trajectory stored from the previous solve (neural-mpc quad_3d_optimizer.py, WebFetch). This is the previous iterate, shifted.
• Approximation order. Paper: first order in the RTI scheme §IV-B. Code takes the order from configuration; the default was not visible in the fetched text (neural-mpc quad_3d_optimizer.py, WebFetch).
• Batched derivatives. Code computes batched Jacobians and Hessians with torch's vmap over jacrev, with and without gradient tracking functorch.py functorch.py. This is the batching described in §IV-A.
• Default activation. The MLP module uses the identity as activation when none is passed ml_perceptron.py. This is not what the paper's networks use, which are not described.

Paper-internal discrepancies, with no code involved:

• The text says the runtime is "below 3ms" for the best networks §VII-A, but the table lists 3.2 ms for N-5-128-Ours tab. III.
• The text says naive optimisation takes "up to 36ms" §VII-A. The table lists 178.7 ms for N-5-128-Naive tab. III. 35.9 ms is N-4-64-Naive.
• Table I gives "up to 11M" parameters for this method tab. I. Table II lists 13M and 12M tab. II, and the text mentions over 150 million parameters at 100 Hz on an RTX 3000 §V fig. 3, which is not in Table II.

Pitfalls

• Section numbers follow the IEEEtran roman numerals used in the compiled paper (e.g. §IV-B, §VII-A). Equation, table and figure numbers were counted from the LaTeX source, in order of appearance, and were not checked against the compiled PDF. Equations 2 and 3 are tagged in an align* with a manual counter increment, so check their numbers against the PDF before citing them externally.
• Speed claim depends on the platform. The 4000x figure is a GPU result on the embedded platform §V fig. 3. The CPU figure in the same paragraph is about 60x. Naive implementations have no GPU row in Table II.
• The approximation is first-order in the RTI scheme. The second-order term is present in the formulation §IV-A eq. 4, but the timing results are for the first-order variant. Do not assume the second-order term costs nothing.
• Error of the approximation is claimed to be second order in step-size terms for the first-order case §IV-A. The paper supports this only with experiments on agile platforms with small δt; no bound is given.
• Unstated inputs for re-implementation: δt, horizon T, cost weights, input bounds, the mapping from rotor thrusts to collective thrust and torque, network activation, weight initialisation, input and output normalisation, epochs and early-stopping patience. Take them from the code only with the caveat in Paper vs code.
• Real-time window is an extra assumption. The 50 Hz and 100 Hz windows are targets the authors chose; the real-world control rate is not stated fig. 3 §VII-B.
• Several network sizes were not run on the real robot. Only N-3-32-Ours and N-3-32-Naive appear in flight tab. IV. The naive version crashed.
• Real-world sample size is small. Five runs per experiment tab. IV; the ground-effect result is from one evaluation with 8 flyovers §VII-B.
• The RTI description is prose. There is no algorithm environment in the paper §IV-B fig. 2. The pseudocode on this card is ours.
• Successor framework. The deprecated ml-casadi repository is what this card cites. For new work the README points to L4CasADi, which includes a real-time mode. Check that project's documentation for the current interface.
• Batched derivatives and BatchNorm. The batched Jacobian path does not support BatchNorm README.md.
• Unlicensed text in the experiment repo. The experiment repo is GPL-3.0. Reimplementing its behaviour from this card is fine; copying code from it is not.

Lineage

• Builds on: the Gaussian-process residual models of Kabzan et al. and Torrente et al. (the latter's quadrotor setup is used for the simplified simulator) §I §VII-A. Neural residual models of Bauersfeld et al. and Saviolo et al. §II. RTI scheme of Diehl et al. §IV-B.
• Compared against: a naive integration of neural dynamics into CasADi, as used by Spielberg et al., Chee et al. and Saviolo et al. §IV-C §VII-A; GP models with 20, 50 and 100 support points tab. III; the nominal model.
• Superseded by: the successor framework L4CasADi, as stated in the ml-casadi README (external, no card).