Paper deep dive
Optimal Control with Learned Critics under Unmodeled State Dependencies
Philipp Schoch, Markus Ryll
Intelligence
Status: not_run | Model: - | Prompt: - | Confidence: 0%
Entities (0)
Relation Signals (0)
No relation signals yet.
Cypher Suggestions (0)
No Cypher suggestions yet.
Abstract
Abstract:Model Predictive Control (MPC) provides a structured and constraint-aware mechanism for decision-making, but its reliance on optimization-friendly analytical dynamics models limits its use in tasks with contacts and other hard-to-model state dependencies. Model-free reinforcement learning avoids explicit modeling assumptions but typically requires large amounts of interaction data. We present a learning-based MPC framework that combines the data efficiency and structure of local model-based planning with learned components that compensate for incomplete dynamics and finite-horizon myopia. The method augments a nominal analytical model with a residual dynamics network that learns missing state-dependent effects from data and combines the resulting planner with a learned action-value critic that injects long-horizon MDP structure into the local iLQR optimization. To make this practical at reinforcement-learning scale, we develop a GPU-accelerated batched iLQR solver that evaluates learned dynamics and critic networks inside the optimal-control loop and solves thousands of trajectory-optimization problems in parallel. The complete system is integrated into a robotics simulator, enabling scalable model-based reinforcement learning under incomplete dynamics. Experiments on biased and incompletely modeled control tasks show that the approach improves closed-loop control performance while preserving the model-based structure needed for efficient constrained trajectory optimization.
Tags
Links
- Source: https://arxiv.org/abs/2610.05359v1
- Canonical: https://arxiv.org/abs/2610.05359v1
Trouble viewing inline? Open PDF directly →
Full Text
71,947 characters extracted from source content.
Expand or collapse full text
Optimal Control with Learned Critics under Unmodeled State Dependencies Philipp Schoch Markus Ryll Affiliation: Autonomous Aerial Systems Lab Affiliation: Technical University of Munich Email: philipp.schoch,markus.ryll@tum.de Abstract Model Predictive Control (MPC) provides a structured and constraint-aware mechanism for decision-making, but its reliance on optimization-friendly analytical dynamics models limits its use in tasks with contacts and other hard-to-model state dependencies. Model-free reinforcement learning avoids explicit modeling assumptions but typically requires large amounts of interaction data. We present a learning-based MPC framework that combines the data efficiency and structure of local model-based planning with learned components that compensate for incomplete dynamics and finite-horizon myopia. The method augments a nominal analytical model with a residual dynamics network that learns missing state-dependent effects from data and combines the resulting planner with a learned action-value critic that injects long-horizon MDP structure into the local iLQR optimization. To make this practical at reinforcement-learning scale, we develop a GPU-accelerated batched iLQR solver that evaluates learned dynamics and critic networks inside the optimal-control loop and solves thousands of trajectory-optimization problems in parallel. The complete system is integrated into a robotics simulator, enabling scalable model-based reinforcement learning under incomplete dynamics. Experiments on biased and incompletely modeled control tasks show that the approach improves closed-loop control performance while preserving the model-based structure needed for efficient constrained trajectory optimization. Keywords: Model Predictive Control, Model-Based Reinforcement Learning, iLQR, Value Function Approximation 1 Introduction Computing optimal control policies for agents acting in stochastic environments is a central challenge in sequential decision making. Formally, this amounts to solving a Markov Decision Process (MDP) [17]: finding a policy that maximizes expected cumulative reward over potentially infinite horizons. The field has produced a rich landscape of approaches, broadly organized into trajectory optimization, model-free reinforcement learning (RL), and model-based reinforcement learning methods [13]. Model Predictive Control (MPC) performs trajectory optimization in a receding-horizon fashion, enabling precise control and explicit satisfaction of state and input constraints. However, MPC requires an optimization-friendly analytical dynamics model, often difficult or impossible to derive for systems involving contacts or other highly nonlinear effects, and breaks down if relevant state variables remain unmodeled. Contact-rich tasks are a representative case, because making and breaking contact introduces discontinuous phase boundaries over which gradient information can become uninformative [16], while classical multi-phase or mixed-integer formulations lead to large nonlinear programs that are brittle and rarely compatible with real-time replanning [7]. A second limitation is the finite horizon itself. MPC applies a finite-horizon surrogate to an infinite-horizon MDP, so a horizon that is too short can make the controller greedy, trap it in local minima, or prevent task completion when intermediate subgoals must be reached before the final objective becomes attainable [6]. A third limitation is that standard MPC formulations operate on explicit state variables and analytical models, making it nontrivial to use high-dimensional sensor observations such as rgb or depth images directly inside the planner. Model-free RL sidesteps these limitations: learned policies are cheap at inference and can consume high-dimensional observations. However, model-free methods are notoriously sample-inefficient, scale poorly with state-space dimensionality, and struggle with long-horizon tasks where reward-sparsity impedes exploration [17, 1]. Model-based RL promises to combine the data efficiency of planning with the flexibility of learning [13]. A compelling paradigm in this space is to use MPC as a critic: the planner provides a structured, locally accurate estimate of the action-value function, which can be used to bootstrap a global Q-approximator which after training supplies long-horizon guidance back to the MPC, allowing even to correct for mismatches in the dynamics model. [3] demonstrates this in the setting of gradient-free sampling-based control. A key remaining challenge is that MPC critics still require an accurate dynamics model, a requirement that fails whenever relevant state variables lack a tractable analytical description, e.g. object contact dynamics in manipulation. In this work, we develop a learning-augmented constrained MPC framework that embeds two learned components directly into the planner. A residual dynamics model augments the nominal analytical model and exposes otherwise unmodeled state dependencies to optimization, while a learned critic injects global MDP information into the finite-horizon objective to reduce myopia and improve guidance under multi-phase cost-structure; unlike in sampling-based MPC [3], the critic enters the Riccati recursion through its first- and second-order derivatives. Following the DDP formulation, we realize this approach with a GPU-accelerated batched ALM-iLQR solver that integrates both learned components inside the optimal-control loop and scales to thousands of parallel environments, enabling efficient model-based reinforcement learning. Core contributions: (i) a constrained MPC formulation that combines residual dynamics learning with a learned critic, and (i) a learning-based DDP-style solver that makes this formulation practical for model-based reinforcement-learning, scaling to thousands of parallel environments. 2 Preliminaries We consider a discrete-time MDP ℳ=(,,l,P,γ,p0)M=(S,A,l,P,γ,p_0) [17], where S is the state space, A the action space, l(x,u)l(x,u) a stage cost, P a stochastic transition model with xk+1∼P(⋅∣xk,uk)x_k+1 P(· x_k,u_k), γ∈(0,1]γ∈(0,1] a discount factor, and p0p_0 the initial state distribution. Let ℙℳπ(⋅∣x,u)P^π_M(· x,u) denote the trajectory distribution induced by (x0,u0)=(x,u)(x_0,u_0)=(x,u) and following π thereafter. The action-value and value functions are Qℳπ(x,u)=ℙπℳ(⋅∣x,u)[∑i=0∞γil(xi,ui)]Q^π_M(x,u)=E_P^π_M(· x,u)\! [ _i=0^∞γ^il(x_i,u_i) ] and Vℳπ(x)=u∼π(⋅∣x)[Qℳπ(x,u)]V^π_M(x)=E_u π(· x)\! [Q^π_M(x,u) ]. The control objective is π∗=argminπx∼p0[Vℳπ(x)]π^*= _π\,E_x p_0\! [V^π_M(x) ]. For any horizon H≥2H≥ 2, the action-value function admits the finite-horizon decomposition Qℳπ(x,u)=ℙπℳ(⋅∣x,u)[∑i=0H−2γil(xi,ui)+γH−1Qℳπ(xH−1,uH−1)].Q^π_M(x,u)=E_P^π_M(· x,u)\! [ _i=0^H-2γ^il(x_i,u_i)+γ^H-1Q^π_M(x_H-1,u_H-1) ]. (1) Classical MPC closes this recursion with a terminal cost, Qℳπ(xH−1,uH−1)≈lT(xH−1,uH−1)Q^π_M(x_H-1,u_H-1)≈ l_T(x_H-1,u_H-1) [8]. While sufficient when the terminal cost is well designed, this local approximation may miss the global structure of the MDP, especially when the objective requires reaching an intermediate subgoal [6]. In the following we work with deterministic planned trajectories to simplify notation; the stochastic setting is retained by taking expectations over the induced trajectory distribution. The paper is structured as follows: Section 3 introduces the blended-critic objective and explains how learning of the residual dynamics function and the learned critic integrates into the framework. Section 4 describes the iLQR solver and the specific realization of the blended-critic objective. Section 5 evaluates the method on a set of simulated control tasks, and Sections 6 and 7 discuss the method’s shortcomings and give an outlook on potential extensions. 3 Optimal Control Objective 3.1 Blending the Learned Critic As established in Eq. (1), MPC estimates QℳπQ^π_M by rolling out a local trajectory. However, the accuracy of this estimate depends heavily on the length of the rollout and well-designed terminal cost function. A global learned critic Q^θ(x,u) Q_θ(x,u), parameterized by a neural network with weights θ, can in principle capture the full MDP structure, but suffers from the usual approximation errors of neural networks, making it both hard to train and hard to integrate into a planner (cf. Section 3.2, cf. [17]). Following [3], we therefore blend both estimates along the H-step planned trajectory (xi,ui)i=0H−1(x_i,u_i)_i=0^H-1. As can be seen later, this will allow us to trade off bias in the nominal model with approximation bias in the learned critic, while reducing variance in the value targets during training. Following the general Bellman recursion Qi=li+γQi+1Q_i=l_i+γ Q_i+1, for any λ∈[0,1]λ∈[0,1], we define the blended Q-estimate recursively for i=0,…,H−2i=0,…,H-2 as Qλ(xi,ui)=(1−λ)Q^(xi,ui)+λ(l(xi,ui)+γQλ(xi+1,ui+1)),Q^λ(x_i,u_i)=(1-λ) Q(x_i,u_i)+λ (l(x_i,u_i)+γ\,Q^λ(x_i+1,u_i+1) ), (2) but contrary to [3], we terminate the recursion at the last step by Qλ(xH−1,uH−1)=(1−λ)Q^(xH−1,uH−1)+λl(xH−1,uH−1)Q^λ(x_H-1,u_H-1)=(1-λ) Q(x_H-1,u_H-1)+λ\,l(x_H-1,u_H-1), as a weighted average between the finite horizon approximation of MPC and the learned critic. Following Appendix A we can rollout this recursion into the explicit form Q0λ=Q^0+λ(∑k=0H−2(λγ)kBk+(λγ)H−1(lH−1−Q^H−1)),Q_0^λ= Q_0+λ ( _k=0^H-2(λγ)^kB_k+(λγ)^H-1 (l_H-1- Q_H-1 ) ), (3) where Bi=li+γQ^i+1−Q^iB_i=l_i+γ Q_i+1- Q_i is the one-step Bellman residual of the learned critic. Setting λ=0λ=0 recovers Q Q directly; λ=1λ=1 recovers the standard MPC rollout; intermediate values let Q Q act as a global shaping term while preserving local model accuracy. In Section 5, we additionally vary the blend weight of the critic at the end of the horizon separately, denoted λf _f; in all other experiments, λf=λ _f=λ. 3.2 Learning the Critic The critic Q^θ Q_θ is trained by stochastic gradient descent on the mean squared Bellman error [17] ℒQ(θ)=(xk,uk,yk)∼[(Q^θ(xk,uk)−yk)2].L_Q(θ)=E_(x_k,u_k,y_k) [ ( Q_θ(x_k,u_k)-y_k )^2 ]. (4) The target yky_k is constructed either from short rollout fragments or from the MPC continuation value Vλ(xk+1)=minuQλ(xk+1,u)V^λ(x_k+1)= _uQ^λ(x_k+1,u) which can be obtained for free after each MPC solve. TD(n) targets. A first option is to use standard n-step targets built from contiguous rollout fragments [17]. This is simple and efficient, and larger n can stabilize the targets through the usual bias-variance trade-off. However, the method becomes sensitive to stale replay data when buffers are small or policy updates are infrequent. MPC-based targets. A second option is to use this continuation value and form the one-step target ykMPC=l(xk,uk)+γVλ(xk+1)y_k^MPC=l(x_k,u_k)+γ V^λ(x_k+1), with the usual truncation at terminal transitions. In the offline variant, the planner stores the successor blended-Q-estimate during data collection and training reuses this cached value [3]. This is fast and benefits from the full planner lookahead, but the targets become stale as the critic changes. In the online variant, the successor value is recomputed during training by solving a fresh MPC problem for each minibatch element. This removes staleness and supports fully off-policy replay, but is typically prohibitively expensive: with iLQR, each batch element requires a converged optimal control solve, and shuffled minibatches prevent amortizing convergence across time steps. Value ensembles. Both target constructions remain bootstrapped and therefore inherit the usual risk of amplifying value-estimation errors [20]. In our setting this is especially problematic because MPC will immediately exploit critic mismatches induced by epistemic uncertainty. To reduce this effect, we follow the spirit of double Q-learning and ensemble bootstrapping [20, 9] and use an ensemble critic with N value heads. During each SGD step, one head is randomly selected for updating, while the remaining N−1N-1 heads provide the target values. Their outputs can be aggregated by a (soft/hard) maximum, minimum, or mean operator; under our cost convention, the maximum is conservative, the minimum is optimistic, and we use the mean as a balanced compromise [9]. At inference time, the ensemble disagreement and the disagreement of the ensembles gradients also provide an uncertainty estimate of the learned critic. This gives us a way to downweigh the λ-contribution of in the MPC in regions where the Value-heads disagree. Contrary to [6], where a softmax (softmin under our cost-convention) is used to encourage exploration, we compute the standard deviation across the normalized value-head outputs (and their gradients) and attenuate the learned-critic contribution by an exponential factor α=e−κσα=e^-κσ where κ is a tuning parameter. 3.3 Learning the Dynamics Residual Many control tasks involve state variables that are directly observable, but for which no sufficiently accurate analytical dynamics model is available. In manipulation, for example, one may observe the pose and velocity of a grasped object without modeling the underlying contact dynamics explicitly. To allow the planner to reason over such quantities, we augment the state with components whose dynamics are incomplete and learn only the missing part from data. Formally, we partition the state as x=[xp⊤,xs⊤]⊤x=[x_p ,x_s ] , where xp∈ℝnpx_p ^n_p denotes primary states with known dynamics and xs∈ℝnsx_s ^n_s denotes secondary states whose evolution must be learned. We then use the composite dynamics x˙=ftot(x,u)=fp(x,u)+fϕ(x,u), x=f_tot(x,u)=f_p(x,u)+I\,f_φ(x,u), (5) where =diag(m1,…,mn)I=diag(m_1,…,m_n) with mi∈0,1m_i∈\0,1\ is a binary mask selecting which state dimensions receive a learned residual correction, while dimensions with trusted dynamics, such as trivial integrators, are left unchanged. We fit the residual from rollout data. Given samples (xk,uk,xk+1)(x_k,u_k,x_k+1), a one-step target is obtained from the finite difference x˙k≈(xk+1−xk)/Δt x_k≈(x_k+1-x_k)/ t; for states containing quaternion components, standard finite differences are replaced by logarithmic-map differences on SO(3)SO(3) so that the targets are well defined despite the q≡−q≡-q sign ambiguity, as described in Appendix H. Because one-step fitting alone does not control the error accumulated over repeated numerical integration inside the planner, we optimize both over mini-batches of the local regression loss and a rollout-consistency loss on contiguous trajectory fragments, ℒ(ϕ)=∑k‖(ftot(xk,uk)−x˙k)‖2,ℒroll(ϕ)=1R∑t=1R‖(xk+tpred−xk+ttrue)/Δt‖2,L(φ)= _k \|I (f_tot(x_k,u_k)- x_k ) \|^2, _roll(φ)= 1R _t=1^R \|I(x_k+t^pred-x_k+t^true)/ t \|^2, (6) where xk+tpredx_k+t^pred is obtained by recursively integrating ftotf_tot from the previous predicted state. This term penalizes compounding rollout error on the residual-enabled state components, rather than only minimizing local one-step errors. The resulting model allows MPC to plan directly in the augmented state space. Training then alternates between collecting rollouts, fitting the residual dynamics, and updating the learned critic. 3.4 Guiding the Critic Integration of Auxiliary Information. Rather than learning directly from raw (x,u)(x,u), both learned modules in the solver have the option to tap into task-specific feature pre-processors, x˙=fp(x,u)+fϕ(gf(x,u,o)),Q^θ(x,u∣o)=Q^θ(gQ(x,u,o)), x=f_p(x,u)+I\,f_φ (g_f(x,u,o) ), Q_θ(x,u o)= Q_θ (g_Q(x,u,o) ), (7) where o denotes optional task context, such as a goal, object pose, or a latent visual embedding z=eψ(o)z=e_ψ(o) produced by a shared visual encoder. The maps gfg_f and gQg_Q may differ and let the learned components operate in task coordinates rather than embodiment-specific raw states. For manipulation tasks, for instance, it can be convenient to compute inside gQg_Q forward-kinematics-derived features such as pcube−pEEp_cube-p_E, pgoal−pcubep_goal-p_cube, and vEEv_E instead of raw joint variables. 4 OCP Solver 4.1 iLQR Approaching the MPC from the MDP perspective (cf. Section 2), the DDP [8] based iLQR approach [19] is a natural choice for our optimal control solver. Given an initial control sequence (ui)i=0H−1(u_i)_i=0^H-1, a forward rollout under xi+1=f(xi,ui)x_i+1=f(x_i,u_i) produces a reference trajectory (xi,ui)i=0H−1(x_i,u_i)_i=0^H-1 and its rollout value V. Following DDP, the backward pass uses the Bellman recursion Qi=li+γQi+1Q_i=l_i+γ Q_i+1 and forms a local quadratic model of the action-value function at each time step, Q(δx,δu)≈12[1δxδu]T[0QxTQuTQxQxxQxuQuQuxQuu][1δxδu].Q(δ x,δ u)≈ 12 bmatrix1\\ δ x\\ δ u bmatrix^\!T bmatrix0&Q_x^T&Q_u^T\\ Q_x&Q_x&Q_xu\\ Q_u&Q_ux&Q_u bmatrix bmatrix1\\ δ x\\ δ u bmatrix. (8) Minimizing this local model yields the affine feedback update δui=ki+Kiδxi,ki=−Quu,i−1Qu,i,Ki=−Quu,i−1Qux,i.δ u_i=k_i+K_iδ x_i, k_i=-Q_u,i^-1Q_u,i, K_i=-Q_u,i^-1Q_ux,i. (9) iLQR drops the second-order dynamics tensor terms of full DDP, yielding the more economic Gauss-Newton approximation [19]. In the backward sweep, we compute the feed-forward terms kik_i, feedback gains KiK_i, and the expected value-change terms used by the line search. The stage derivatives, dynamics Jacobians, and Hessians can be batch-computed before this sequential sweep, which allows us to GPU-parallelize the heavy computations related to auto-diff gradient and Hessian evaluation beforehand and keep the backward-sweep loop as lightweight as possible. We regularize QuuQ_u using the scheme of Tassa et al. [18], with an aggressive schedule for decreasing or increasing the regularizer depending on line-search success. The forward line search rolls out candidate controls u^i=ui+αki+Kidx(x^i,xi) u_i=u_i+α k_i+K_i\,dx( x_i,x_i) from x^0=x0 x_0=x_0, where α is the line-search step size [18]. The state displacement dx(⋅,⋅)dx(·,·) is Euclidean for vector states and uses the manifold distance for quaternion states, as described in Appendix H. We accept a candidate when the actual rollout improvement agrees with the expected local-model improvement, using either an Armijo condition or the ratio rule of Tassa et al. [18]; otherwise α is reduced and the rollout is repeated. 4.2 Augmented Lagrangian Method To handle equality and inequality constraints, we wrap the iLQR primal solve in an Augmented Lagrangian Method (ALM). Let τ=(x0,u0,…,xH−1,uH−1)τ=(x_0,u_0,…,x_H-1,u_H-1) denote the planned trajectory. We write the constrained finite-horizon problem with a generic additive objective J(τ)=∑i=0H−1l~i(xi,ui)J(τ)= _i=0^H-1 l_i(x_i,u_i), where l~i l_i will later be instantiated by the blended MPC stage objective: minτ _τ J(τ)=∑i=0H−1l~i(xi,ui) J(τ)= _i=0^H-1 l_i(x_i,u_i) (10) subject to to x0=x, x_0=x, xi+1=f(xi,ui), x_i+1=f(x_i,u_i), gi(xi,ui)≤0, g_i(x_i,u_i)≤ 0, hi(xi,ui)=0,i=0,…,H−1. h_i(x_i,u_i)=0, i=0,…,H-1. Let ci=[hi⊤,gi⊤]⊤c_i=[h_i ,g_i ] collect constraint residuals and let μi=[μh,i⊤,μg,i⊤]⊤ _i=[ _h,i , _g,i ] collect the corresponding multipliers. For fixed multipliers and penalties, the iLQR subproblem uses the augmented Lagrangian stage objective [12, 4], LiAL(xi,ui,μi,ρi)=l~i(xi,ui)+μi⊤ci(xi,ui)+12ci(xi,ui)⊤Iρ,ici(xi,ui).L_i^AL(x_i,u_i, _i, _i)= l_i(x_i,u_i)+ _i c_i(x_i,u_i)+ 12c_i(x_i,u_i) I_ρ,ic_i(x_i,u_i). (11) The diagonal penalty matrix is Iρ,i=diag(ρh,i,mi⊙ρg,i)I_ρ,i=diag( _h,i,m_i _g,i), where mij=0m_ij=0 when gij(xi,ui)<0g_ij(x_i,u_i)<0 and μg,ij=0 _g,ij=0, and mij=1m_ij=1 otherwise. Thus, inactive inequality constraints do not contribute residuals, gradients, or Hessians to the local model. ALM alternates between a primal and a dual step. In the primal step, iLQR minimizes ∑i=0H−1LiAL _i=0^H-1L_i^AL with μi _i and Iρ,iI_ρ,i fixed, using the rollout and line-search procedure of Section 4.1. After the inner solve returns the primal trajectory (xi⋆,ui⋆)(x_i ,u_i ), the dual step updates μh,i+=μh,i+ρh,ihi(xi⋆,ui⋆),μg,i+=max(0,μg,i+ρg,igi(xi⋆,ui⋆)). _h,i^+= _h,i+ _h,ih_i(x_i ,u_i ), _g,i^+= (0, _g,i+ _g,ig_i(x_i ,u_i ) ). (12) For constraints that remain violated, the penalty is increased as ρi+=min(Φρi,ρmax) _i^+= ( _i, _max) with Φ>1 >1. When the multiplier update vanishes, the residuals approach zero, and stationarity of the augmented subproblem reduces to the KKT stationarity conditions of the original constrained problem [10]. 4.3 Blended-iLQR The constrained blended-MPC objective is obtained by substituting the ALM stage objective LiAL(xi,ui,μi,ρi)L_i^AL(x_i,u_i, _i, _i) into the blended recursion of Section 3.1. Define the marginalized continuation value Vi+1λ(xi+1)=minui+1Qiλ(xi+1,ui+1)V_i+1^λ(x_i+1)= _u_i+1Q_i^λ(x_i+1,u_i+1). For fixed ALM multipliers and penalties, the constrained recursion is Qiλ(xi,ui) Q_i^λ(x_i,u_i) =(1−λ)Q^(xi,ui)+λ(LiAL(xi,ui,μi,ρi)+γVi+1λ(xi+1)) =(1-λ) Q(x_i,u_i)+λ (L_i^AL(x_i,u_i, _i, _i)+γ V_i+1^λ(x_i+1) ) (13) =λLiAL(xi,ui,μi,ρi)+(1−λ)Q^(xi,ui)⏟L~iλ(xi,ui)+λγVi+1λ(f(xi,ui)). = λ L_i^AL(x_i,u_i, _i, _i)+(1-λ) Q(x_i,u_i)_ L_i^λ(x_i,u_i)+λγ V_i+1^λ(f(x_i,u_i)). At the last control step, the continuation is omitted and QH−1λ=L~H−1λQ_H-1^λ= L_H-1^λ. Thus the blended AL-iLQR has the same dynamic-programming form as standard iLQR, but with blended ALM stage-objective L~iλ L_i^λ and effective discount λγλγ. The backward pass therefore uses the same quadratic expansion as Eq. (8); only the derivatives of the local objective change. In particular, the constrained stage derivatives combine the ALM derivatives with the learned critic derivatives, e.g., L~x,iλ=λ(LiAL)x+(1−λ)Q^x,L~u,iλ=λ(LiAL)u+(1−λ)Q^u, L_x,i^λ=λ(L_i^AL)_x+(1-λ) Q_x, L_u,i^λ=λ(L_i^AL)_u+(1-λ) Q_u, (14) and analogously for L~xx,iλ L_x,i^λ, L~uu,iλ L_u,i^λ, and L~ux,iλ L_ux,i^λ. The full expansions are given in Appendix B. The remaining recursion is unchanged, except that the next-state value contribution is scaled by λγλγ. 5 Results We evaluate the proposed framework in four complementary experiments. We begin with a runtime benchmark to verify our solver is amenable to the large-throughput requirement of reinforcement learning. We test the solver on randomly generated optimal-control problems to characterize how our implementation scales with batch size, hardware, and the method of derivative computation. Next, we isolate the effect of the critic blend on a point-mass task, where the optimal critic is available in closed form. We then study a quadrotor task with a biased dynamics model and an unmodeled suspended payload, which tests whether joint residual-dynamics learning and critic learning can compensate for model bias, capture missing state couplings, and remain compatible with quaternion kinematics. Finally, we consider Franka cube stacking as a contact-rich manipulation task to assess whether the framework can plan over incompletely modeled contact interactions and benefit from visual and haptic information for grasp detection. Figure 1: Parallel environments used for training and evaluation. Left: The stack-cube task with a 6-dof manipulator. The robot must pick up the target cube (blue) and stack it on top of the goal cube (red). Right: The quadrotor task, where the drone must hover at a target position while carrying an unmodeled, suspended payload. Both tasks are simulated using Genesis [2]. Speed Benchmarking. The experimental setup and detailed timing curves are given in Appendix C. The main result is that GPU execution becomes decisive once learned dynamics and critic models are part of the solve: on the tested hardware, batched GPU execution is roughly one to two orders of magnitude faster than CPU execution over the relevant batch-size range. In the purely non-learning setting, CPU execution is competitive only for the smallest batches, after which GPU execution overtakes it as batching exposes more parallel work. Figure 4 further shows that the amortized runtime per problem decreases strongly up to about 10410^4 parallel QPs before saturating, and that automatic differentiation remains close to analytical derivatives in this benchmark. Isolating the Effect of the Critic Blend. To separate the effect of the blend from critic approximation errors, a force-controlled 2-D point mass must reach a fixed goal position and come to rest. For this linear-quadratic task, the optimal infinite-horizon discounted Q∗Q^* and V∗V^* follow in closed form from the discrete algebraic Riccati equation and take the place of the learned critic (Appendix D). We sweep the blend for horizons H∈2,5,10,20H∈\2,5,10,20\ in two settings: blending Q∗Q^* at all stages with λf=λ _f=λ, and blending only the terminal value with λf _f, while λ=1λ=1 at all other stages. Figure 2 shows the results. Without the critic (λ=λf=1λ= _f=1), MPC overshoots the goal for H∈2,5,10H∈\2,5,10\ and fails the task, whereas H=20H=20 succeeds. Lowering λf _f towards zero recovers the success rate and cost of the optimal policy at all horizons. Blending at all stages reaches the same performance already at values of λ closer to 11. We interpret this as the stage-wise critic terms supplying long-range information along the whole planned trajectory, in addition to the terminal value. These results rely on the exact Q∗Q^*; Section 6 discusses why learned critics restrict λ to values close to 11. Figure 2: Point-mass reach task with the analytical LQR critic Q∗,V∗Q^*,V^*. Columns: success rate and mean episode cost. Top: all-stage blend, λ=λfλ= _f swept jointly. Bottom: terminal-only blend, λ=1λ=1 and λf _f swept. Overcoming Biased Dynamics and Unmodeled State-Dependencies. In this experiment, a quadrotor starts in a random pose and must hover to a target position while carrying an unmodeled payload suspended by a rod. With all learning components disabled, the nominal MPC baseline reaches a success rate of only 0.050.05 and a mean cost of 2.6e42.6e4, with most failures caused by crashes. The learned critic variant, which combines residual-dynamics learning with the blended critic objective, raises the success rate to 0.790.79 and lowers the mean cost to 4.7e34.7e3. As an adaptive-control baseline, an idealized incremental nonlinear dynamic inversion (INDI) controller [15] without sensor noise or delay reaches a higher success rate of 0.940.94, but a higher mean cost of 5.6e35.6e3, since it overshoots more and settles more slowly around the goal (Appendix E.1). In our setting, learning the dynamics model alone is not sufficient to solve the task reliably. We attribute this to the short look-ahead of 11s for the MPC, not capturing the fatal swing-up motion of the pendulum early enough, and the approximation bias of the learned dynamics model, requiring careful alignment of planning frequency, integration scheme, and look-ahead distance to trade off approximation bias with real-time-ness. However, we observe that combining the learned dynamics model with the learned critic, which injects infinite-horizon MDP information into the MPC, yields more conservative approaches towards the goal, keeping the oscillations bounded and avoiding crashing. The learned critic can drive the system away from high action-value regions that lead to failure, including the ones with high pendulum deflection. This observation is broadly in line with the remarkable success of model-free policy-gradient methods, which achieve high task performance not by meticulously modeling the underlying dynamics but by directly learning mitigation strategies that avoid task failure [11]. To further demonstrate the sample-efficiency claim, we establish a PPO baseline [14] with the hyperparameters of Table 6 reaching comparable asymptotic performance but requiring training for 1.5 hours on 1e41e4 parallel environments, collecting roughly 1e81e8 environment transitions (collected at 200 Hz, the policy frequency), whereas our learned-critic achieves the reported performance after roughly 30 minutes using 10510^5 transitions (collected at 20 Hz, the planner frequency, cf. Appendix G). Figure 3 shows how this improvement emerges during training: as λ is decreased, the learned critic gradually takes over a larger share of the objective, correcting the nominal behavior and improving both success rate and cost. Detailed task and training settings are given in Appendix E.1, whereas empirical evidence for the efficacy of our approach and further insights into the learned critic are given in Appendix E.2. Figure 3: Quadrotor learning progress while decreasing the critic blend parameter λ. Early updates behave like nominal MPC and fail under the biased incomplete model (cf. Section 3.1); once the learned critic is allowed to shape the objective more strongly, the success rate increases and the cost drops. Learning to Plan over Contacts Using Visual and Haptic Information. In this experiment, a 6-DOF manipulator must pick up a cube and stack it precisely on top of another cube using only a simple kinematic robot model and no explicit contact model. Unsurprisingly, a naive grasp-uninformed nominal MPC with a single global cost function fails completely (success rate of 00, stuck in pre-grasp pose) since it cannot reason over the grasp. Following Appendix F.1, we add to the state vector the cube’s position and a grasp-score, which is initialized from an external grasp detection using tactile sensors, and replace the global cost with a multi-phase cost function. This modification raises the nominal baseline to 0.3610.361 success rate. Further learning a residual to the grasp-score state lifts the success rate to moderate 0.4010.401. Employing a learned critic, which is conditioned on visual depth and haptic observations (cf. Figure 9) improves task performance significantly to 0.940.94 success rate. We observe more delicate and better-informed recovery behaviors than with nominal MPC, including corrections after failed grasp attempts and poor approaches that tip over the target cube. More insights into the task setup and further discussions of the results are given in Appendix F.1 and Appendix F.2. 6 Limitations A first limitation is that the solver relies on first- and second-order derivatives of the learned residual and critic models, and optimization remains sensitive to irregular derivatives. This restricts the blend: with the exact critic of Section 5, lower values of λ consistently recover the performance of the optimal policy, whereas with learned critics only modest reductions from λ=1λ=1, roughly in the range 0.980.98 to 0.90.9, outperformed the nominal baselines in our experiments. Since the weight of the nominal rollout decays as (γλ)k(γλ)^k along the horizon (Eq. (3)), low λ removes the model information that would compensate imperfect critic derivatives, particularly at long horizons (cf. Appendix F.2). We therefore start training at λ=1λ=1, lower λ only once the critic has acquired sufficient state-space coverage, and mitigate irregular derivatives by disagreement-based attenuation of the critic, smoothing, and clipping. Nevertheless, poor gradients can still steer the optimizer toward uninformative regions of the state space, degrade data collection early-on, and ultimately stall training; stronger training-time regularization, efficient off-policy target construction and priority-sampling from the replay-buffer should therefore be investigated. A second limitation concerns real-time deployment on a single system. While all solver computations are batched across environments, and costs and derivatives are evaluated in batch over all time steps, the forward rollout, the Riccati recursion, and the line search remain sequential over the horizon, which precludes replanning at rates above 100Hz100\,Hz in the current implementation. However, iLQR returns the time-indexed local feedback policy u^i=ui+Kidx(x^i,xi) u_i=u_i+K_i\,dx( x_i,x_i) around the planned trajectory, so the planner can run at a low rate while this policy closes the loop at the control rate. 7 Conclusion To conclude, this paper presented a learning-augmented constrained MPC framework for planning under model mismatch. While dynamics learning enabled the planner to account for directly observable but incompletely modeled state variables, the main driver of the observed performance gains was the integration of a learned critic into the gradient-based optimization framework. By providing corrective descent directions during the backward recursion, the critic injected long-horizon behavioral information that is inaccessible to a finite-horizon MPC objective alone. In this way, the method mitigated the inherent myopia of MPC by leveraging value information from the underlying infinite-horizon MDP. Across the reported runtime and control experiments, the framework scaled to thousands of parallel optimal control problems and improved performance under biased and incomplete dynamics models. Future work will focus on addressing the aforementioned limitations and on stronger critic pre-training, including representation learning from video demonstrations. Appendix A Bellman-Residual Form of the Blended Objective We have H control steps, hence the planned rollout contains the state-action pairs (xi,ui)i=0H−1(x_i,u_i)_i=0^H-1 and terminal state xHx_H. For compactness, write Qi=Q(xi,ui)Q_i=Q(x_i,u_i), Q^i=Q^(xi,ui) Q_i= Q(x_i,u_i), and li=l(xi,ui)l_i=l(x_i,u_i). For i=0,…,H−2i=0,…,H-2, the blended recursion in the main text reads Qiλ=(1−λ)Q^i+λ(li+γQi+1λ),Q_i^λ=(1-λ) Q_i+λ (l_i+γ Q_i+1^λ ), where Q^i Q_i is the learned Q-estimate and QiλQ_i^λ is the MPC-blended Q estimate at timestep i. At the last control step we close the recursion with QH−1λ=(1−λ)Q^H−1+λlH−1.Q_H-1^λ=(1-λ) Q_H-1+λ l_H-1. To expose the relation to temporal-difference style corrections, define the deviation of the blended objective from the learned critic by δi=Qiλ−Q^i. _i=Q_i^λ- Q_i. Subtracting Q^i Q_i from both sides of the recursion gives δi=(1−λ)Q^i+λ(li+γQi+1λ)−Q^i=λ(li+γQi+1λ−Q^i). _i=(1-λ) Q_i+λ (l_i+γ Q_i+1^λ )- Q_i=λ (l_i+γ Q_i+1^λ- Q_i ). Now add and subtract γQ^i+1γ Q_i+1 inside the parentheses: δi=λ(li+γQ^i+1−Q^i+γ(Qi+1λ−Q^i+1)). _i=λ (l_i+γ Q_i+1- Q_i+γ(Q_i+1^λ- Q_i+1) ). This separates the local Bellman residual of the learned critic from the downstream blended correction. Define the one-step Bellman residual Bi=li+γQ^i+1−Q^i.B_i=l_i+γ Q_i+1- Q_i. Using δi+1=Qi+1λ−Q^i+1 _i+1=Q_i+1^λ- Q_i+1, we obtain the recursion δi=λ(Bi+γδi+1), _i=λ(B_i+γ _i+1), At the last control step, the terminal approximation gives δH−1=QH−1λ−Q^H−1=λ(lH−1−Q^H−1). _H-1=Q_H-1^λ- Q_H-1=λ (l_H-1- Q_H-1 ). For any intermediate step we get δi=λ(∑k=0H−i−2(λγ)kBi+k+(λγ)H−i−1δH−1). _i=λ ( _k=0^H-i-2(λγ)^kB_i+k+(λγ)^H-i-1 _H-1 ). Setting i=0i=0 and substituting back Q0λ=Q^0+δ0Q_0^λ= Q_0+ _0 yields Q0λ=Q^0+λ(∑k=0H−2(λγ)kBk+(λγ)H−1(lH−1−Q^H−1)).Q_0^λ= Q_0+λ ( _k=0^H-2(λγ)^kB_k+(λγ)^H-1 (l_H-1- Q_H-1 ) ). Hence, the blended objective can be interpreted as the learned critic plus a geometrically weighted accumulation of local Bellman-residual corrections along the planned trajectory. In particular, when the learned critic satisfies the Bellman relation well along the planned rollout, the residual terms BiB_i are small and the blended objective stays close to Q^0 Q_0; when the local model-based rollout disagrees with the learned critic, the residual terms accumulate and correct the objective accordingly. Appendix B Quadratic Q Expansion for blended ALM-iLQR For the ALM-iLQR the gradients and hessians for the quadratic Q expansion in each timestep read: Qx Q_x =L~x+fx⊤Vx′, = L_x+f_x V_x , Qu Q_u =L~u+fu⊤Vu′, = L_u+f_u V_u , Qxx Q_x =L~xx+fx⊤Vxx′fx+Vx′⋅fxx, = L_x+f_x V_x f_x+V_x · f_x, Quu Q_u =L~uu+fu⊤Vuu′fu+Vx′⋅fuu, = L_u+f_u V_u f_u+V_x · f_u, Qux Q_ux =L~ux+fu⊤Vxx′fx+Vx′⋅fux. = L_ux+f_u V_x f_x+V_x · f_ux. with the learned Q-function terms, and using a locally-linear approximation of the constraint function ci≈ci0+cx,iδxi+cu,iδuic_i≈ c_i^0+c_x,iδ x_i+c_u,iδ u_i, L~x L_x =λ(LiAL)x+(1−λ)Q^x, =λ(L_i^AL)_x+(1-λ) Q_x, (LiAL)x (L_i^AL)_x =lx+cx,i⊤(μi+Iρ,ici), =l_x+c_x,i ( _i+I_ρ,ic_i), L~u L_u =λ(LiAL)u+(1−λ)Q^u, =λ(L_i^AL)_u+(1-λ) Q_u, (LiAL)u (L_i^AL)_u =lu+cu,i⊤(μi+Iρ,ici), =l_u+c_u,i ( _i+I_ρ,ic_i), L~xx L_x =λ(LiAL)xx+(1−λ)Q^xx, =λ(L_i^AL)_x+(1-λ) Q_x, (LiAL)xx (L_i^AL)_x ≈lxx+cx,i⊤Iρ,icx,i, ≈ l_x+c_x,i I_ρ,ic_x,i, L~uu L_u =λ(LiAL)uu+(1−λ)Q^uu, =λ(L_i^AL)_u+(1-λ) Q_u, (LiAL)uu (L_i^AL)_u ≈luu+cu,i⊤Iρ,icu,i, ≈ l_u+c_u,i I_ρ,ic_u,i, L~ux L_ux =λ(LiAL)ux+(1−λ)Q^ux, =λ(L_i^AL)_ux+(1-λ) Q_ux, (LiAL)ux (L_i^AL)_ux ≈lux+cu,i⊤Iρ,icx,i. ≈ l_ux+c_u,i I_ρ,ic_x,i. the next-state value-function gradient Vx′V_x and hessian Vxx′V_x are carried over in the recursion. Appendix C Runtime Scaling Benchmark To isolate solver throughput, we benchmark the batched iLQR implementation on random control instances with an 8-dimensional state, a 4-dimensional action, a planning horizon of 32 steps, and a planning time step of 0.05s0.05\,s. The learned dynamics model and learned critic are fully connected MLPs with three hidden layers of width 256 and ELU activations. Each timing run measures a single iLQR iteration, comprising one forward rollout, one backward pass, up to six regularization trials, and up to eight line-search evaluations. We discard two warmup runs and report averages over 20 measured runs. The benchmark is executed on an NVIDIA RTX 5080 mobile GPU (175 W) and an AMD Ryzen 9 9955HX CPU with 64 GB RAM. Figure 4: Runtime scaling of one batched iLQR iteration with analytical and autodiff derivatives on CPU and GPU. GPU execution becomes decisive once learned models are part of the solve, while autodiff remains close to analytical derivatives; runtime per QP decreases with batching up to about 10410^4 parallel problems before flattening on the tested hardware. The figure separates total solve time from amortized runtime per QP. For the learned-model setting, the GPU curves stay well below the CPU curves across the relevant batch sizes, yielding roughly one to two orders of magnitude lower runtime on the tested hardware. For the analytical setting, CPU execution is competitive only at the smallest batches, after which GPU execution overtakes it as batching exposes more parallel work. The amortized runtime decreases strongly up to about 10410^4 parallel QPs and then saturates, which in our setup coincides with the practical hardware limit of the tested mobile GPU. Appendix D Point-Mass Experiment Details The point mass has mass m=2kgm=2\,kg, state x=[p⊤,v⊤]⊤∈ℝ4x=[p ,v ] ^4 with planar position p and velocity v, and is driven by a planar force u∈[−20,20]2Nu∈[-20,20]^2\,N with p˙=v p=v and v˙=u/m v=u/m. The stage cost is ℓ(x,u)=βp∥p−p⋆∥22+βv∥v∥22+βu∥u∥22 (x,u)= _p p-p _2^2+ _v v _2^2+ _u u _2^2 with βp=1 _p=1, βv=0.05 _v=0.05, and βu=0.1 _u=0.1. Linearizing the discretized dynamics at Δtplan=0.1s t_plan=0.1\,s gives (A,B)(A,B), and the discounted algebraic Riccati equation with γ=0.95γ=0.95 gives P, so that V∗(x)=e⊤PeV^*(x)=e Pe and Q∗(x,u)=ℓ(x,u)+γV∗(Ae+Bu)Q^*(x,u)= (x,u)+γ\,V^*(Ae+Bu) in the error coordinates e=x−[p⋆⊤,0]⊤e=x-[p ,0] . An episode of 10s10\,s counts as successful once the position error falls below 0.1m0.1\,m and the speed below 0.05m/s0.05\,m/s. Each sweep runs 10310^3 parallel environments and holds each blend value for one episode, λ∈1,0.99,0.98,0.97,0.95,0.93,0.9,0.8,…,0λ∈\1,0.99,0.98,0.97,0.95,0.93,0.9,0.8,…,0\ for the all-stage sweep and λf∈1,0.9,…,0 _f∈\1,0.9,…,0\ for the terminal-only sweep. Appendix E Quadrotor Experiment Details E.1 Task Setup Dynamics Model. The nominal MPC model uses the rigid-body state and control given by x x =[pWBWqWBvWBWωWB],u=[azBαWB], = bmatrixp_WB^W\\ q_WB\\ v_WB^W\\ _WB^B bmatrix, u= bmatrixa_z^B\\ _WB^B bmatrix, The virtual controls are subject to componentwise box constraints, umin≤u≤umax,azB∈[−20,40]ms−2,αWB∈[−40,40]3rads−2,u_ ≤ u≤ u_ , a_z^B∈[-20,40]\,m\,s^-2, _WB^B∈[-40,40]^3\,rad\,s^-2, equivalently written in the inequality-constraint form as hmin(u)=umin−u≤0h_ (u)=u_ -u≤ 0 and hmax(u)=u−umax≤0h_ (u)=u-u_ ≤ 0. For the quaternion block, the solver uses the reduced local coordinates and projection matrix xred=[pWBWϕWBvWBWωWB],E=[I0000G(qWB)0000I0000I],x_red= bmatrixp_WB^W\\ _WB\\ v_WB^W\\ _WB^B bmatrix, E= bmatrixI&0&0&0\\ 0&G(q_WB)&0&0\\ 0&0&I&0\\ 0&0&0&I bmatrix, so the quaternion state is optimized in a three-dimensional tangent coordinate ϕWB∈ℝ3 _WB ^3 rather than directly in qWB∈S3q_WB∈ S^3. The dynamics are x˙=f(x,u)=[vWBW12qWB⊗[0ωWB]gW+R(qWB)[00azB]αWB]. x=f(x,u)= bmatrixv_WB^W\\ 12\,q_WB bmatrix0\\ _WB^B bmatrix\\ g^W+R(q_WB) bmatrix0\\ 0\\ a_z^B bmatrix\\ _WB^B bmatrix. where gW=[00−g]g^W= bmatrix0\\ 0\\ -g bmatrix and R(qWB)R(q_WB) maps body-frame accelerations into the world frame. To inject structured model mismatch, the simulator additionally exposes a swing-load state with roll- and pitch pendulum coordinates, xaug=[xθ˙],θ∈ℝ2,x^aug= bmatrixx\\ θ\\ θ bmatrix, θ ^2, corresponding to a 0.172kg0.172\,kg suspended load (0.15kg0.15\,kg payload, 0.02kg0.02\,kg rod, 0.002kg0.002\,kg bearing) attached to a 0.15m0.15\,m rod. These state variables are observed during training but are not modeled in the nominal rigid-body dynamics above; their effect on the translational and rotational motion must therefore be captured from data alone. The model bias m¯=m/1.5 m=m/1.5 and J¯B=1.5JB J_B=1.5J_B is introduced in the the actuator model v˙WBW v_WB^W =gW+m¯−1R(qWB)fB(ω2), =g^W+ m^-1R(q_WB)f_B(ω^2), ω˙WB ω_WB^B =J¯B−1(τB(ω2)−ωWB×J¯BωWB), = J_B^-1\! ( _B(ω^2)- _WB^B× J_B _WB^B ), which computes the the squared rotor-speed commands ωi2 _i^2 from the virtual controls u using the wrench model fB(ω2) f_B(ω^2) =[00kf∑j=14ωj2],τB(ω2)=∑j=14rj×[00kfωj2]+[00kd∑j=14sjωj2], = [ smallmatrix0\\ 0\\ k_f _j=1^4 _j^2 smallmatrix ], _B(ω^2)= _j=1^4r_j× [ smallmatrix0\\ 0\\ k_f _j^2 smallmatrix ]+ [ smallmatrix0\\ 0\\ k_d _j=1^4s_j _j^2 smallmatrix ], with kfk_f and kdk_d the force and drag coefficients, rjr_j the rotor-arm vectors, and sj∈−1,+1s_j∈\-1,+1\ the propeller spin directions. Hence the benchmark combines parameter bias on the modeled rigid-body states with additional unmodeled state dependencies from the swing-load coordinates. Cost Function The quadrotor stage cost is a quadratic hover-tracking objective, where (⋅)⋆(·) denotes the hover goal state, ℓ(x,u)=12(CLOSE (x,u)= 12 ( βp∥pWBW−pWBW,⋆∥22+βq∥rotvec(qWB⋆⊗qWB−1)∥22 _p p_WB^W-p_WB^W, _2^2+ _q (q_WB q_WB^-1) _2^2 +βv∥vWBW−vWBW,⋆∥22+βω∥ωWB−ωWB,⋆∥22 + _v v_WB^W-v_WB^W, _2^2+ _ω _WB^B- _WB^B, _2^2 OPEN+βu∥u∥22), + _u u _2^2 ), where βp _p, βq _q, βv _v, βω _ω, and βu _u are positive cost weights. INDI Baseline. INDI [15] computes control increments from measured accelerations and thereby compensates unmodeled forces and moments, including those of the suspended payload, without an explicit model of them. We evaluate an idealized variant, in which these measurements are free of noise and delay, on the same task and with the same cost function as the MPC variants. E.2 Experimental Insights Failure mode of nominal MPC. Figure 5 compares representative closed-loop trajectories of nominal MPC and the blended-critic controller under the biased incomplete model (random seeds). The nominal controller initially reduces the position error, but it does so with the wrong inertial model and without explicit payload dynamics. This injects excessive rotational motion into the suspended load, which amplifies and destabilizes the quadrotor. Around the end of the episode, the nominal controller tips over, loses altitude and crashes. The learned critic avoids this failure: by learning to favor a lower-velocity approach, it keeps the pendulum motion bounded, avoiding catastrophic tip-over even under presence of severe cross-couplings. Figure 5: Two closed-loop trajectories under the biased incomplete quadrotor model. Solid lines denote the trajectories stabilized by the learned critic controller; dashed lines denote the nominal MPC. Different random initial states are shown for the two controllers. The nominal controller initially moves toward the goal but excites the suspended load, which then lead to tip-over late into the flight. The learned controller approaches the goal more conservatively, keeping the pendulum motion bounded and stabilizing the system throughout the episode. Monte Carlo analysis. The single-rollout picture is confirmed by the Monte Carlo envelopes in Figure 6, computed from N=300N=300 randomized closed-loop runs. Under the learned controller, the mean trajectories converge toward hover and the state dispersion contracts after the initial transient. Under nominal MPC, the dispersion grows late in the episode as the quadrotor tips over and crashes during breaking. The nominal mean also develops a clear downward drift in z. This is the aggregate version of the same failure mode seen in Figure 5. Blended Critic MPC Nominal MPC Figure 6: Monte Carlo envelopes over N=300N=300 randomized rollouts. Left: blended-critic MPC. Right: nominal MPC. The learned controller keeps the mean close to the hover target and prevents late-episode variance growth. The nominal controller develops large dispersion and a systematic loss of altitude, consistent with repeated crash trajectories. Local structure of the learned critic. Figures 7 and 8 show local x-y slices of the learned critic. Each slice is generated from one state-action sample taken during MPC optimization: only the quadrotor position coordinates pxp_x and pyp_y are swept, while the remaining state and action entries are fixed. The plots therefore describe the value density around this specific rollout condition; they are not global projections of the learned Q function, which is fitted over the full 2121-dimensional state-action input. Figure 7 shows that the four value heads converge to almost identical spatial value distributions in this slice. The critic is therefore certain about the correction signal used at this rollout state. Figure 8 compares this learned correction to the nominal stage cost. The critic remains consistent with the stage cost in the important failure case: it does not create spurious low-value basins far away from the goal and remains locally smooth. At the same time, the value minimum is shifted relative to the instantaneous stage-cost minimum. For the fixed velocity, attitude, and payload state in this snapshot, positions closer to the goal have higher predicted continuation cost; the same non-position states would be preferable farther back along the approach direction. This is the correction to account for the increased momentum the nominal model cannot express. The critic surface is smooth over the plotted region, which is required for stable iLQR derivatives. The final panel shows the uncertainty-based scaling derived from value-head variance, which is used to weight the critic against the nominal MPC objective (cf. Section 3.2). Figure 7: Local x-y slices of the four learned value heads at one MPC optimization snapshot. All state and action entries except planar position are held fixed. The heads agree on the spatial value distribution, indicating low ensemble uncertainty for the critic correction at this rollout state. Figure 8: Local critic analysis for the same slice as Figure 7. The learned value remains consistent with the stage cost by avoiding low-value regions far from the goal, but shifts the local optimum according to the fixed velocity, attitude, and payload state. The critic is smooth in the plotted region, and the value-head variance defines the guidance scale used to blend the critic with the nominal MPC objective. Appendix F Cube Stacking Experiment Details F.1 Task Setup Dynamics Model. The manipulation benchmark asks the robot to pick up a target cube and place it on top of a support cube. The planner uses the reduced state and control x x =[qarmwgpWTWg],u=[q˙armw˙g],x˙=f(x,u)=[q˙armw˙gvWEW(qarm,q˙arm)0], = bmatrixq_arm\\ w_g\\ p_WT^W\\ g bmatrix, u= bmatrix q_arm\\ w_g bmatrix, x=f(x,u)= bmatrix q_arm\\ w_g\\ g\,v_WE^W(q_arm, q_arm)\\ 0 bmatrix, where qarmq_arm are arm joints, wgw_g is gripper width, pWTWp_WT^W is target-cube position, and g∈[0,1]g∈[0,1] is a learned continuous grasp gate. Thus the nominal model transports the cube kinematically only in proportion to the current grasp estimate: for g≈0g≈ 0 the target cube is effectively decoupled from the end effector, whereas for g≈1g≈ 1 it follows the commanded end-effector motion. The end-effector translational velocity is computed from differentiable forward kinematics of the manipulator. Let pWEW=FKE(qarm)p_WE^W=FK_E(q_arm) denote the end-effector position obtained from the robot’s batched forward-kinematics model. The corresponding geometric Jacobian is JEW(qarm)=∂pWEW∂qarm,vWEW(qarm,q˙arm)=JEW(qarm)q˙arm.J_E^W(q_arm)= ∂ p_WE^W∂ q_arm, v_WE^W(q_arm, q_arm)=J_E^W(q_arm) q_arm. In the implementation, the forward kinematics and geometric Jacobians are evaluated in batch across all environments and are fully differentiable. Contact dynamics are not modeled explicitly; instead, the missing interaction effects must be absorbed by the learned residual and the learned value function. Cost Function. Let pWAWp_WA^W denote the support-cube center and define the desired stacked target position pWGW=pWAW+[0,0,1.2sz]⊤p_WG^W=p_WA^W+[0,0,1.2s_z] , with sxs_x and szs_z the cube width and height. We implement a multi-phase cost function, where the pregrasp-term and the transport term are weighted by the grasping gate ℓ(x,u)= (x,u)= βp∥pWTW−pWEW∥2+βpclip(1.8sx−pWT,zW,0,1.8sx) _p p_WT^W-p_WE^W _2+ _p\,clip(1.8s_x-p_WT,z^W,0,1.8s_x) +gβp∥pWGW−pWEW∥2+(1−g)sg[βp∥pWGW−pWEW∥2] +g\, _p p_WG^W-p_WE^W _2+(1-g)\,sg\! [ _p p_WG^W-p_WE^W _2 ] −βg(ρ(pWTW)w˙g−(1−ρ(pWTW))gw˙g)+βu∥q˙arm∥22+βgw˙g2, - _g (ρ(p_WT^W)\, w_g- (1-ρ(p_WT^W) )g\, w_g )+ _u q_arm _2^2+ _g w_g^2, where βp _p, βg _g, and βu _u are positive cost weights, ρ(pWTW)=exp(−∥pWGW−pWTW∥2/(0.4sx))ρ(p_WT^W)= \! (- p_WG^W-p_WT^W _2/(0.4s_x) ), and sg[⋅]sg[·] denotes the stop-gradient evaluation used in the implementation. Preprocessed Inputs. The residual dynamics model receives the feature vector ϕf=[x⊤u⊤(fLW)⊤(fRW)⊤gzimg⊤]⊤, _f= bmatrixx &u &(f_L^W) &(f_R^W) &g&z_img bmatrix , that is, the full state, control, left and right finger contact forces, the scalar grasp score, and optional image features zimgz_img. The value model receives ϕQ=[(pWGW)⊤x⊤u⊤(fLW)⊤(fRW)⊤zimg⊤]⊤, _Q= bmatrix(p_WG^W) &x &u &(f_L^W) &(f_R^W) &z_img bmatrix , so it additionally gets the goal location but does not take the grasp score as a separate input in the current implementation. F.2 Experimental Insights The ablations in Tables 1, 2, and 3 support the findings summarized in Section 5. On this task, the main benefit comes from critic learning with informative observations, whereas residual dynamics learning yields only moderate improvement. Effect of planning horizon. Table 1 shows that nominal MPC works best with a medium length optimization horizon. Very short horizons do not result in a sufficiently strong grasp, keeping the robot stuck in the approach phase, while very long horizons degrade due to unsteadiness in the iLQR backward pass when optimizing over the grasp. The learned critic improves success across the full horizon range and already reaches 0.6850.685 success rate at H=5H=5, where nominal MPC barely succeeds. This supports the main-text claim that the critic gives iLQR consistently a descent direction aligned with the long-horizon stacking objective, leading to more directed behavior and more robust grasping. The learned critic also degrades at long horizons; unlike the exact critic of Section 5, it supplies imperfect gradients and curvature at every stage, while the weight of the nominal rollout decays as (γλ)k(γλ)^k, and both effects can destabilize the backward pass over long horizons (cf. Section 6). Table 1: Influence of iLQR horizon length on cube-stacking performance at Δtctrl=0.05s t_ctrl=0.05\,s. Nominal MPC requires an intermediate horizon, whereas the learned critic maintains strong performance even for short horizons. H Nominal SR Nominal cost Critic SR Critic cost 2 0.00 – 0.03 278 5 0.116 280 0.685 266 10 0.318 254 0.885 223 20 0.361 251 0.941 220 30 0.298 253 0.018 335 50 0.0 411 0.02 421 Effect of visual and haptic inputs. Table 2 shows a consistent gain when depth and haptic observations are added to the critic input. Using only state, action, and goal features yields a success rate of 0.230.23. Adding depth raises this to 0.5430.543, and adding haptic information further raises it to 0.9410.941. Figure 9 illustrates the transport phase together with the predicted cube and end-effector trajectories available to the controller. This matches the main-text result: visual and haptic observations stabilize the critic by reducing ambiguity about grasp quality and task phase, which is the dominant source of uncertainty in this benchmark. Figure 9: Cube-stacking phases. (1) Approach and Grasp; (2) Pickup and Transport; (3) Align and Drop; Depth (B) and haptic observations (contact-probes in gripper fingers visualized by red markers, figure A1..3) improve grasp estimation. The predicted cube trajectory (C2, blue) matches the predicted end-effector trajectory (C2, green) during the transport phase (2). Table 2: Influence of critic input features on cube-stacking performance. Adding depth and haptic observations improves success rate monotonically, indicating that richer observations improve grasp-phase inference and critic stability. Critic input SR Mean cost [x,u,pgoal][x,u,p_goal] 0.230 238.0 [x,u,pgoal,depth][x,u,p_goal,depth] 0.543 234.0 [x,u,pgoal,depth,haptic][x,u,p_goal,depth,haptic] 0.941 257.0 Effect of residual dynamics learning. Table 3 shows only minor benefit from learning residual dynamics on this task. The learned variants remain close to or below the baseline success rate and consistently increase the mean cost. This supports the hypothesis stated in the results section: modeling the exact dynamics of the cube is not the bottleneck in this task. Once the task phase is inferred robustly, the phased nominal model is already sufficient to plan approach, lift, and transport, whereas additional residual dynamics introduces optimization non-convexity that hurts iLQR. In practice, the added residual model occasionally induced erratic motions, such as brief gripper opening or slight end-effector jerks. During the transport phase, when the end effector had to maintain a firm grasp, such events often led to dropping the object, reducing success rate and increasing cost. Table 3: Influence of residual dynamics learning on cube-stacking performance. Residual variants do not improve significantly over the nominal baseline and often increase cost. Experiment SR Mean cost Baseline 0.394 234.7 g˙←[x,u,contact] g←[x,u,contact] 0.341 275.6 g˙←[contact,depth] g←[contact,depth] 0.401 278.2 g˙←[x,u,contact,depth] g←[x,u,contact,depth] 0.323 273.4 x˙←[x,u,contact] x←[x,u,contact] 0.361 296.4 Appendix G Solver and Learning Settings Tables 4, 5, and 6 summarize the shared solver and learning configuration, the task-specific settings used in the quadrotor and cube-stacking experiments, and the PPO hyperparameters. Table 4: Shared solver, optimization, and data-collection settings used in both tasks. Setting Value Planning time step Δtplan t_plan 0.05s0.05\,s Horizon H 2020 Parallel environments 100100 iLQR iterations 88 Maximum line-search trials 99 Maximum regularization updates 1010 Executed control steps per replan 33 Residual-model initialization Gaussian, standard deviation 10−310^-3 Residual rollout loss 1010 steps, weight 0.10.1, midpoint integration Critic initialization Gaussian, standard deviation 10010^0 Value-target construction TD(5) Variance-based λ attenuation κ=4.0κ=4.0 Blend schedule 1.0→0.91.0→ 0.9 Update interval every 2020 control steps Transitions per update 2×1032× 10^3 Mini-batches per update 1010 Critic replay buffer 10410^4 Dynamics replay buffer 10710^7 Table 5: Task-specific timing and network settings. Setting Quadrotor Cube stacking Episode length 5s5\,s 8s8\,s Control time step Δtctrl t_ctrl 0.005s0.005\,s 0.01s0.01\,s Residual model ELU MLP, 33 layers, width 128128 ELU MLP, 33 layers, width 256256 Critic trunk ELU MLP, 33 layers, width 128128 ELU MLP, 33 layers, width 256256 Value ensemble 44 heads, 11 layer, width 128128, mean aggregation 44 heads, 11 layer, width 256256, mean aggregation Table 6: PPO hyperparameters. Setting Value Parallel Environments 10410^4 Learning epochs 55 Mini-batches 44 Clipping parameter 0.20.2 Discount factor γ 0.990.99 GAE parameter λ 0.950.95 Value-loss coefficient 1.01.0 Entropy coefficient 0.010.01 Learning rate 10−310^-3 Desired KL 0.010.01 Appendix H Quaternion Projection Often in aerospace and robotics, one is interested in optimizing over singularity-free rotation parametrizations. To make the solver accessible to a broad range of tasks, we require native support for state vectors that mix Euclidean states with quaternion components. During the forward rollout, quaternion states are renormalized after numerical integration, and during line search, the displacement dx(x^i,xi)dx( x_i,x_i) uses Euclidean differences on vector blocks and local rotation errors on quaternion blocks. In the backward pass, quaternion derivatives are projected to local tangent coordinates because direct perturbations in the four-dimensional quaternion coordinates introduce radial components that violate ‖q‖=1\|q\|=1 and make the local quadratic model inconsistent after renormalization [5]. We therefore map derivatives from the quaternion representation q to a local angle-axis coordinate ϕ∈ℝ3φ ^3 using the geometric Jacobian G(q)G(q). For a quaternion-dependent function h, this gives hxred=E⊤hxh_x^red=E h_x and hxxred≈E⊤hxxEh_x^red≈ E h_xE, where the Hessian projection omits second-order curvature terms in the spirit of the iLQR Gauss-Newton approximation. For dynamics, input and output states lie in different tangent spaces, so fx,ired=Ei+1⊤fx,iEif_x,i^red=E_i+1 f_x,iE_i and fu,ired=Ei+1⊤fu,if_u,i^red=E_i+1 f_u,i [5]. Appendix E.1 gives a concrete quadrotor example, and in the following we will detail the construction of G(q)G(q). Given a current unit quaternion q and a local Rodrigues perturbation ψ∈ℝ3ψ ^3, define δq(ψ)=11+‖ψ‖2[1ψ],q+(ψ)=q⊗δq(ψ).δ q(ψ)= 1 1+\|ψ\|^2 bmatrix1\\ ψ bmatrix, q^+(ψ)=q δ q(ψ). For a function h:S3→ℝph:S^3 ^p, the local rotational Taylor expansion is h(q+(ψ))=h(q)+∂h(q⊗δq(ψ))∂ψ|ψ=0ψ+(‖ψ‖2).h(q^+(ψ))=h(q)+ ∂ h(q δ q(ψ))∂ψ |_ψ=0ψ+O(\|ψ\|^2). Using left quaternion multiplication q⊗δq=L(q)δq δ q=L(q)δ q, the projected Jacobian is ∂h∂ψ=∂h∂qL(q)∂δq∂ψ|ψ=0. ∂ h∂ψ= ∂ h∂ qL(q) ∂δ q∂ψ |_ψ=0. For the Cayley parametrization, ∂δq∂ψ|ψ=0=[01×3I3]and thereforeG(q)=L(q)[01×3I3], ∂δ q∂ψ |_ψ=0= bmatrix0_1× 3\\ I_3 bmatrix therefore G(q)=L(q) bmatrix0_1× 3\\ I_3 bmatrix, (15) which gives the result in [5]. with q=[w,v⊤]⊤q=[w,v ] and quaternion left-multiplication matrix L(q)=[w−v⊤vwI3+[v]×].L(q)= bmatrixw&-v \\ v&wI_3+[v]_× bmatrix. For the Hessian, we reuse the first-order projection matrix from Eq. (15), yielding the projected Hessian (without 2nd-order curvature) hϕϕ=G(q)⊤hqq(q)G(q).h_φ=G(q) h_q(q)G(q). References [1] A. Aubret, L. Matignon, and S. Hassas (2019) A survey on intrinsic motivation in reinforcement learning. ArXiv abs/1908.06976. External Links: Link Cited by: §1. [2] G. Authors (2024) Genesis: a generative and universal physics engine for robotics and beyond. External Links: Link Cited by: Figure 1. [3] M. Bhardwaj, S. Choudhury, and B. Boots (2020) Blending MPC & value function approximation for efficient reinforcement learning. CoRR abs/2012.05909. External Links: Link, 2012.05909 Cited by: §1, §1, §3.1, §3.1, §3.2. [4] T. A. Howell, B. E. Jackson, and Z. Manchester (2019) ALTRO: a fast solver for constrained trajectory optimization. In 2019 IEEE/RSJ International Conference on Intelligent Robots and Systems (IROS), Vol. , p. 7674–7679. External Links: Document Cited by: §4.2. [5] B. E. Jackson, K. Tracy, and Z. Manchester (2021) Planning with attitude. IEEE Robotics and Automation Letters 6 (3), p. 5658–5664. External Links: Document Cited by: Appendix H, Appendix H. [6] K. Lowrey, A. Rajeswaran, S. M. Kakade, E. Todorov, and I. Mordatch (2018) Plan online, learn offline: efficient learning and exploration via model-based control. CoRR abs/1811.01848. External Links: Link, 1811.01848 Cited by: §1, §2, §3.2. [7] T. Marcucci, J. Umenberger, P. A. Parrilo, and R. Tedrake (2021) Shortest paths in graphs of convex sets. CoRR abs/2101.11565. External Links: Link, 2101.11565 Cited by: §1. [8] D. MAYNE (1966) A second-order gradient method for determining optimal trajectories of non-linear discrete-time systems. International Journal of Control 3 (1), p. 85–95. External Links: Document, Link, https://doi.org/10.1080/00207176608921369 Cited by: §2, §4.1. [9] O. Peer, C. Tessler, N. Merlis, and R. Meir (2021) Ensemble bootstrapping for q-learning. CoRR abs/2103.00445. External Links: Link, 2103.00445 Cited by: §3.2. [10] (2006) Penalty and augmented lagrangian methods. In Numerical Optimization, p. 497–528. External Links: ISBN 978-0-387-40065-5, Document, Link Cited by: §4.2. [11] J. Peters and J. A. Bagnell (2010) Policy gradient methods. In Encyclopedia of Machine Learning, C. Sammut and G. I. Webb (Eds.), p. 774–776. External Links: ISBN 978-0-387-30164-8, Document, Link Cited by: §5. [12] M. J. D. Powell (1969) A method for nonlinear constraints in minimization problems. External Links: Link Cited by: §4.2. [13] R. Reiter, J. Hoffmann, D. Reinhardt, F. Messerer, K. Baumgärtner, S. Sawant, J. Boedecker, M. Diehl, and S. Gros (2025) Synthesis of model predictive control and reinforcement learning: survey and classification. External Links: 2502.02133, Link Cited by: §1, §1. [14] C. Schwarke, M. Mittal, N. Rudin, D. Hoeller, and M. Hutter (2025) RSL-rl: a learning library for robotics research. arXiv preprint arXiv:2509.10771. Cited by: §5. [15] E. J. J. Smeur, Q. Chu, and G. C. H. E. de Croon (2016) Adaptive incremental nonlinear dynamic inversion for attitude control of micro air vehicles. Journal of Guidance, Control, and Dynamics 39 (3), p. 450–461. External Links: Document, Link, https://doi.org/10.2514/1.G001490 Cited by: §E.1, §5. [16] H. J. T. Suh, M. Simchowitz, K. Zhang, and R. Tedrake (2022) Do differentiable simulators give better policy gradients?. CoRR abs/2202.00817. External Links: Link, 2202.00817 Cited by: §1. [17] R. S. Sutton and A. G. Barto (2018) Reinforcement learning: an introduction. Second edition, The MIT Press. External Links: Link Cited by: §1, §1, §2, §3.1, §3.2, §3.2. [18] Y. Tassa, T. Erez, and E. Todorov (2012) Synthesis and stabilization of complex behaviors through online trajectory optimization. In 2012 IEEE/RSJ International Conference on Intelligent Robots and Systems, Vol. , p. 4906–4913. External Links: Document Cited by: §4.1, §4.1. [19] E. Todorov and W. Li (2005) A generalized iterative lqg method for locally-optimal feedback control of constrained nonlinear stochastic systems. Proceedings of the 2005, American Control Conference, 2005., p. 300–306 vol. 1. External Links: Link Cited by: §4.1, §4.1. [20] H. van Hasselt, A. Guez, and D. Silver (2015) Deep reinforcement learning with double q-learning. CoRR abs/1509.06461. External Links: Link, 1509.06461 Cited by: §3.2.