teach-an-engineer · optimization · embedded control

Inside acados

How an open-source NMPC solver squeezes a nonlinear optimal-control problem into a 1 millisecond budget: discretization, Gauss–Newton SQP, structure-exploiting QP solving, and the real-time iteration.

1Why acados exists: control is a live optimization problem

Motivation first. Model predictive control (MPC) is the "do the right thing by imagining the future" controller: at every sample, it solves a constrained optimal-control problem (OCP) over a receding horizon, applies only the first input, and repeats. On paper this is beautiful — constraints, cost shaping, and model knowledge all enter natively. In practice there is one brutal requirement: the OCP must be solved every sample period, faster than it elapses, on modest hardware, or the loop closes on stale data and the controller is useless.

For a linear system with quadratic cost and no constraints, MPC is one matrix factorization. Everything interesting — road-vehicle path planning, quadrotor flight, robot manipulation, chemical plants — is nonlinear MPC (NMPC): the dynamics and often the constraints are nonlinear, so no closed form exists and you must iterate a numerical solver inside the control loop. A horizon with 50 stages of a 6-state robot yields a nonlinear program (NLP) with thousands of decision variables that general-purpose NLP solvers solve in seconds — three to four orders of magnitude too slow for a 1–10 ms loop.

acados (from a-cados Code Generation, its predecessor; by the syscop group at University of Freiburg) exists to close that gap. Its answer is not one clever trick but a pipeline of five ideas, each of which strips work that a generic solver does unnecessarily:

Figure 1 — the acados thesis: five layers of structure exploitation
Nonlinear OCP over a receding horizon — thousands of variables, must re-solve every 1–10 ms
1 · Multiple-shooting discretization + fast integrators → a structured NLP, not an ODE
2 · Gauss–Newton SQP: one linearization per sample replaces a "solve the NLP"
3 · The QP subproblem is block-tridiagonal → Riccati / condensing (HPIPM), never a dense factorization
4 · Real-time iteration: warm-start from last cycle's shifted plan, stop after one QP
All of it on BLASFEO panel-major linear algebra + codegen — 0.2–10 ms per cycle on embedded CPUs
Generic solvers treat MPC as a black-box NLP. acados treats it as a control problem and refuses to do work that the problem's structure makes unnecessary. Every remaining section expands one box.

Throughout, we carry one concrete example — the same one the demo at the end runs live. A cart-pole (pole balanced on a force-driven cart): state $x=(p,\dot p,\theta,\dot\theta)$ with $n_x=4$, scalar force input $u$, $n_u=1$, one pole constraint absent, target $x^{\rm ref}=(0.8,\,0,\,0,\,0)$ (pole upright, cart at 0.8 m), horizon $N=16$ stages at $\Delta t = 50$ ms. That is already a nonlinear OCP — and acados solves its RTI cycle in well under a millisecond on a laptop core. Hold on to the numbers $n_x=4,\ n_u=1,\ N=16$; they reappear in every section.

Takeaway: acados is fast because it never solves "an NLP" — it repeatedly solves the tiny, structured pieces that a warm-started Newton method on a multiple-shooting MPC problem actually needs.

2The problem acados solves: the OCP-NLP and multiple shooting

Motivation. Before algorithms, we need the exact data structure acados consumes — its entire speed story is "exploit the shape of this problem", so we must see the shape. acados accepts an optimal control problem in Bolza form with a phase structure, discretized into $N+1$ nodes:

$$ \begin{aligned} \min_{\{x_i,u_i\}} \;& \tau_N\,\ell_N(x_N) \;+\; \sum_{i=0}^{N-1}\tau_i\,\ell_i(x_i,u_i)\\[2pt] \text{s.t.}\quad & x_{i+1} = f_i(x_i,u_i), \qquad i=0,\dots,N-1\\ & \underline{z}_i \le h_i(x_i,u_i) \le \bar z_i,\quad i=0,\dots,N-1, \qquad \underline{z}_N \le h_N(x_N)\le \bar z_N\\ & x_0 = \hat x_0 \qquad\text{(measured state enters here)} \end{aligned} $$

Symbols: $x_i\in\mathbb{R}^{n_x}$ state at node $i$; $u_i\in\mathbb{R}^{n_u}$ input held on stage $i$; $f_i$ the integrator image of the ODE $\dot x = f_c(x,u)$ over the stage length $\tau_i$; $\ell_i$ running cost, $\ell_N$ terminal cost; $h_i$ general path constraints (box bounds $\underline x_i\le x_i\le\bar x_i$, $\underline u_i\le u_i\le\bar u_i$ are a special case). Bounds can be made soft with slack variables and a linear penalty — the exact-L1 trick that lets SQP never fail from an infeasible constraint instead of a transient disturbance. In our example, $\ell_i = (x_i-x^{\rm ref})^\top Q (x_i-x^{\rm ref}) + r\,u_i^2$ with $Q=\mathrm{diag}(30,2,200,4)$, $r=0.02$.

Why multiple shooting — the discretization that defines everything else

Discretizing the OCP into an NLP has a classic trade-off, and acados bets hard on one side:

Figure 2 — single vs. multiple shooting: who is a decision variable
single shooting — only the $u_i$ are unknowns; states are chained outputs (sensitivities $\partial x_i/\partial u_j$ multiply along the chain)
start
$x_0=\hat x_0$
$\phi_1(u_0)$
simulated
$\phi_2(u_0,u_1)$
simulated
$\phi_3(u_0,u_1,u_2)$
simulated
$\cdots = x_N$
multiple shooting — every $x_i$ is an independent variable; dashed links are the soft-to-zero constraints $x_{i+1}-f_i(x_i,u_i)=0$, Newton’s job is to close them
$x_0$
measured
$x_1$ ?
$x_2$ ?
$x_3$ ?
$\cdots\; x_N$ ?
Single shooting simulates everything from $x_0$, so the optimizer’s sensitivity $\partial x_i/\partial u_j$ is a product of $i$ Jacobians — exponentially ill-conditioned for unstable dynamics (the pole is one: unstable pole at $\approx 4\ \mathrm{s^{-1}}$) — and the Newton matrix becomes dense. Multiple shooting guesses the whole trajectory up front and lets Newton pull the dashed continuity gaps to zero; each gap touches only three blocks, so the Jacobian stays block-banded — the structure §5 sells for an $O(N)$ factorization.

Concretely, the decision vector acados hands its QP layer is

$$ Z \;=\; (x_0,u_0,\;x_1,u_1,\;\dots,\;x_{N-1},u_{N-1},\;x_N)\ \in\ \mathbb{R}^{(N+1)n_x + N n_u}, \qquad V_i(Z) \;=\; f_i(x_i,u_i) - x_{i+1} \;=\; 0 , $$

where $V_i$ is the shooting-gap (continuity) constraint at stage $i$: "the integrator says the next state is $f_i(x_i,u_i)$; node $x_{i+1}$ agrees only if the gap vanishes." For our cart-pole: $Z \in \mathbb{R}^{17\cdot 4 + 16} = \mathbb{R}^{84}$ — 84 primal variables, 64 equality gaps. Small by NLP standards; still hopeless to solve "exactly" at 20 Hz with a black-box method — but perfectly sized for what §3–§5 do next.

Takeaway: acados turns the OCP into one big NLP whose variables are all states and inputs, tied together by stage-wise gap constraints — trading a chain of simulation sensitivities for a sparse, banded constraint Jacobian.

3Integrators: turning $\dot x = f_c(x,u)$ into the $f_i$ of the NLP

Motivation. Section 2's NLP needs $f_i(x_i,u_i)$ — the state after $\tau_i$ seconds of continuous dynamics — plus its Jacobians, thousands of times per control cycle. That is numerical ODE integration, and the choice of integrator is a genuine algorithmic decision in acados, not plumbing: it sets accuracy, stage length (hence $N$ and QP size), and how much derivative work costs.

acados ships several integrators, organized around a hierarchy of how nonlinear the dynamics are in the input:

Figure 3 — the integrator menu, from explicit to exploit-structure

ERK — explicit Runge–Kutta

$x_{i+1} = x_i + \tau\sum_j b_j k_j$, with stages $k_j = f_c(\cdot)$ evaluated downstream only.
Fixed order 1–8 (Butcher tableaus, code-generated). Sensitivities by forward sensitivity ODEs alongside. Stable only for $\tau\lesssim$ stability limit — fine for fast robots. acados' default.

IRK — implicit RK (Radau IIA)

Stages satisfy $k_j = f_c(\cdot)$ implicitly: each stage is a small Newton solve sharing one LU.
A-stable: accurate at large $\tau$, so $N$ stages can cover a long horizon. Newton iterations inside each $f_i$ — more CPU per stage, far fewer stages. (3-stage Radau IIA: order 5.)

GNSF — generalized nonlinear sensitivity function

Split $\dot x = A x + B u + \Phi_{\mathrm{lin}}(x,u) + \Phi_{\mathrm{nl}}(z)$ with the nonlinear part depending only on the low-dimensional channel $z = C_{\mathrm{nl}}x + D_{\mathrm{nl}}u$, $n_z\ll n_x$.
Integrate the small nonlinear part, then recover state and sensitivity propagation through it with matrix operations ("sensitivity function"). Cost scales with $n_z$, not $n_x$. For dynamics whose nonlinearity lives in a few coordinates: $\sin\theta$-only models, drag terms, etc.

CVODES — adaptive BDF/Adams (external SUNDIALS)

Variable order/step, adjoint sensitivities, root events.
For stiff or stiff-switching models where hard real-time bounds on your integrator are unrealistic. Variable work per call = awkward for certified worst-case timing.

Discrete / external

Provided $f_i(x,u)$ in C, or via symbolic expression graph (CasADi), auto-differentiated once at model-build time.
Exact maps, learned models, or hand-tuned F12-style physics. acados asks only "give me $f_i$ and $\partial f_i$ fast".

Simulation-only mode

Integrate forward with no sensitivities (also for testing/validation of a plant model).
Sensitivities are the expensive half — skip them when you're only simulating.
The dividing question: how nonlinear in $u$, and how stiff? Explicit RK if dynamics are smooth and slow relative to $\tau$; implicit RK if you need long stages; GNSF when the nonlinearity is low-dimensional — then sensitivity cost collapses to a few matrix–matrix products.

Why sensitivities are the real product. SQP (next section) never uses $f_i$ alone — it uses the pair $\bigl(f_i,\ \partial f_i\bigr)$, where $\partial f_i = [\,\partial_x f_i \;\; \partial_u f_i\,] \in \mathbb{R}^{n_x\times(n_x+n_u)}$. For an RK scheme with $s$ stages, forward sensitivity propagation differentiates each Butcher stage analytically:

$$ X_1 = x_i,\quad G_1 = [\,I\;\;0\,]; \qquad X_j = x_i + \tau\textstyle\sum_{l\lt j} a_{jl}k_l,\quad G_j = [\,I\;\;0\,] + \tau\textstyle\sum_{l\lt j} a_{jl}\frac{\partial k_l}{\partial(x_i,u_i)}; \qquad \frac{\partial k_j}{\partial(x_i,u_i)} = D_1 f_c(X_j,u_i)\,G_j + D_2 f_c(X_j,u_i)[\,0\;\;I\,], $$

so the Jacobian rides along "for free" during the same $k$-evaluations — one pass over the stage graph, both outputs ($X_j$ = stage abscissa, $G_j=\partial X_j/\partial(x_i,u_i)$; with input lifting the $[\,0\;\;I\,]$ picks the stage's own input block). Each $D_1 f_c$ evaluation is where code-generated, expression-graph-differentiated C code (from CasADi) earns its keep: it is a straight-line instruction sequence, no expression-tree walking at run time. This is exactly why the skill's whole acados lineage — ACADO → acados — invested in symbolic/algorithmic differentiation: 90% of one SQP linearization is just "evaluate $f_i$ and $\partial f_i$ for all $N$ stages".

Takeaway: an acados integrator sells (image, Jacobian) as one product; GNSF and IRK exist purely to buy that pair more cheaply per unit of accuracy (longer stages) or per unit of nonlinearity (low-dimensional nonlinear split).

4Gauss–Newton SQP: one Newton step per sample is the whole algorithm

Motivation. With the OCP-NLP of §2 in hand, the standard nonlinear-programming answer is sequential quadratic programming (SQP): repeatedly replace the NLP by a QP that matches it to second order at the current guess, solve the QP exactly, move to the QP's solution, repeat. The reason this is the right family for MPC — not interior point on the NLP, not penalty methods — is that the QP subproblem inherits the multi-stage block-tridiagonal structure, which is the only structure any solver can exploit for $O(N)$ work (§5). acados' OCP solvers (the ocp_nlp module: full-step SQP, and the "SQP with 1 iteration" RTI variant) are SQP implementations in this classic sense.

4.1 What the QP subproblem is

At iterate $(\bar x_i,\bar u_i)$, SQP minimizes a quadratic model of the Lagrangian subject to linearized constraints. acados makes one deliberate approximation, and its name tells you the trade: Gauss–Newton (the "G" in many acados papers). Write the stage Lagrangian (multipliers $\lambda_i$ for the gap $V_i$, $\mu_i$ for path constraints):

$$ \mathcal L = \sum_{i=0}^{N-1} \tau_i\ell_i + \tau_N\ell_N + \sum_{i=0}^{N-1} \lambda_i^\top\big(f_i(x_i,u_i)-x_{i+1}\big) + \text{constraint terms}. $$

A full Newton step needs $\nabla^2_{ZZ}\mathcal L$, which drags in second derivatives of the dynamics $\partial^2 f_i$ — expensive (differentiate the codegen again, and evaluate twice the work) and often the numerically worst-behaved term (it can make the QP indefinite). Gauss–Newton simply drops those, keeping the curvature the problem actually needs:

$$ H_i \;=\; \tau_i\,\nabla^2 \ell_i(\bar x_i,\bar u_i) \;+\; \underbrace{\textstyle\sum_{j=1}^{n_x} \lambda_{i,j}\,\nabla^2 f_{i,j}(\bar x_i,\bar u_i)}_{\text{dropped by Gauss--Newton}} \;\;\approx\;\; \tau_i\begin{bmatrix} Q_i & 0\\[2pt] 0 & R_i \end{bmatrix} \;\succeq\; 0, $$

i.e. the exact stage Hessian block is "cost curvature + multiplier-weighted second derivatives of the integrator map $f_{i,j}$"; GN keeps only the cost curvature (for $h_i$, the same logic keeps $\mu_i\,\nabla^2 h_{i,j}$ only if you enable con_hessian_type = EXACT_HESSIAN — acados offers both, defaulting to GN/none).

This is the same approximation Levenberg–Marquardt uses for least squares: when residuals can be driven small (constraints can be satisfied — a soft-constrained MPC plan usually can), GN converges locally at full rate, and you have guaranteed a convex QP ($R_i\succ 0$) that a Riccati recursion can solve. That convexity guarantee is not a nicety; it is what makes the whole pipeline robust on a production cycle.

The QP subproblem — decision variables the increments $\Delta x_i, \Delta u_i$ plus slacks $s$:

$$ \min_{\Delta Z,\,s}\;\tfrac12 \Delta Z^\top H\,\Delta Z + g^\top \Delta Z + \text{slack costs} \quad\text{s.t.}\quad \underbrace{A_{\rm tridiag}}_{\text{linearized gaps }x_{i+1} = A_i x_i + B_i u_i + d_i}\Delta Z = r_{\rm gap},\quad \underbrace{C\,\Delta Z + s}_{\text{linearized }h_i} \in [\underline{\ } ,\overline{\ }], $$

where the affine gap term $d_i = f_i(\bar x_i,\bar u_i) - A_i \bar x_i - B_i \bar u_i$ and the residual $r_{\rm gap}$ carry the current continuity error — this is Newton pulling the dashed arrows of Figure 2 to zero.

Figure 4 — one SQP iteration on the cart-pole NLP
① Linearize all stages
$f_i \to (A_i,B_i,d_i)$, cost $\to(Q_i,R_i)$, constraints $\to$ Jacobians
= integrator sensitivity pass ×$N$
② Solve the OCP-QP
GN Hessian, gaps as equalities
→ $O(N)$ Riccati / condensing
gives $\Delta Z^\star$ and feedback gains $K_i$
③ Line search / update
$(\bar x,\bar u) \mathrel{+}= \alpha\,\Delta Z^\star$
merit function $=$ Lagrangian + $\sigma\|r_{\rm gap}\|_1$
④ Stop? Convergence measure
KKT residuals of QP→NLP: $\|r_{\rm gap}\|$, dual feas., complementarity
classical SQP: iterate ①–④ to tolerance (offline NLP use)
iterate to tolerance …
… or (RTI) exit after ONE loop — §6
Two personalities of the same solver: with the loop iterated to tolerance, acados is a legitimate direct-multiple-shooting trajectory optimizer for offline OCPs. With the loop cut to one pass per sampled measurement — the dotted exit — it becomes an online feedback controller, which is what §6 is about.

4.2 The merit-function subtlety (why SQP ≠ "just Newton")

Newton on the KKT system can diverge on a nonlinear program (the quadratic model lies; a step can be huge). acados therefore takes a line search on a filter or $\ell_1$ merit function — Lagrangian plus a penalty on constraint violation — damping only the step length $\alpha$, never the direction. Because the GN Hessian is uniformly $\succeq$ the true Lagrangian Hessian (it's the cost Hessian; dynamics curvature dropped), directions are descent directions and a step $\alpha\le 1$ always improves the merit — the mechanism behind SQP's famously reliable global behavior on MPC problems.

Takeaway: SQP turns the NLP into a sequence of convex QPs; the Gauss–Newton choice (drop $\partial^2 f$) is what makes every one of those QPs both cheap to build and guaranteed convex — trading global Newton exactness for per-iteration cost and robustness, a trade MPC's warm-start regime easily repays.

5The QP engine: block-tridiagonal structure, Riccati, and condensing

Motivation. The SQP step's QP is where acados spends most of its milliseconds, and it is where the multiple-shooting bet of Fig. 2 pays out. The equality-constrained QP from §4 has, in condensed-in-time variables $(\Delta x_i,\Delta u_i)$, a KKT matrix whose sparsity is block-tridiagonal: stage $i$ talks only to $i-1$ and $i+1$, because the gap $V_i$ couples only $(x_i,u_i,x_{i+1})$. A generic sparse Cholesky finds this too, but with hundreds of microseconds of index bookkeeping and poor cache behavior. acados hands the QP to HPIPM, which knows the multi-stage form by construction and offers two rival ways to exploit it. Understanding both — and the dial between them — is understanding modern embedded MPC.

5.1 The ladder from "spread out" to "condensed"

Figure 5 — one QP, three equivalent shapes, three solvers
OCP-QP (all $x_i,u_i$ free) — KKT is block-tridiagonal
after partial condensing ($N_c$ macro-stages)
after full condensing — dense in $(\Delta x_0,\Delta u_\bullet)$
Riccati / RTPI
eliminate stages backwards, $N$ tiny $n_u\times n_u$ solves
$O\!\big(N\,n_x^3 + N n_x^2 n_u\big)$, never forms a big matrix
Factorized $N_c$ blocks
Riccati over $N_c$ blocks + small dense tail
best of both: parallelizable macro-stages
Dense Cholesky/LDL$^\top$
one $O(n_{\rm cond}^3)$ factorization + box-constrained IPM
best when $n_{\rm cond}=(N{+}1)n_x{+}Nn_u$ stays small
solve directly
condense groups
condense all
Condensing = substituting each gap constraint away: $\Delta x_{i+1} = A_i\Delta x_i + B_i\Delta u_i + d_i$ turns $x_1..x_N$ into affine functions of $(\Delta x_0, u_\bullet)$, shrinking the KKT system at the price of dense blocks (cost $O(N^2 n_x^3)$ build-up).
Riccati = the same substitution, but organized backward and never materializing the dense result: it recurses a cost-to-go $P_i$, so complexity stays linear in $N$. It's the discrete-time LQR Riccati recursion — applied to the QP, gains included for free.
All three solve the identical QP to identical solution — only the work organization changes. The acados knobs are qp_solver ∈ {PARTIAL_CONDENSING_HPIPM, FULL_CONDENSING_HPIPM, FULL_CONDENSING_QPOASES, …} plus qp_solver_cond_N (the number of macro-stages $N_c$); the RTI path in C additionally exposes the Riccati-based primal-dual IPM (legacy option name RTI_RICATI_DP). Riccati: small $n_x$ (robot), many stages. Condensing: tiny $n_x$ (cart-pole, $n_x{=}4$!) or few stages. Middle ground for everything else. Inequality handling differs too: Riccati gets a reduced-space primal-dual IPM over box/soft constraints; condensed forms can also go to an active-set solver (qpOASES) warm-started in the active set from the last cycle.

5.2 The Riccati recursion, concretely

Because you will see these exact recursions in hpipm's ocp_qp_ipm and in HPIPM's docs — and they are three lines of logic — here they are on the unconstrained core of the QP (constraints handled by the IPM wrapper that repeatedly feeds Riccati modified-cost data). Deviation coordinates $\delta x_i$, stage cost $(\delta x_i+e_i)^\top Q_i(\delta x_i+e_i) + \delta u_i^\top R_i\delta u_i$, linearized dynamics $\delta x_{i+1} = A_i\delta x_i + B_i\delta u_i + d_i$:

$$ \underbrace{K_i = (R_i + B_i^\top P_{i+1}B_i)^{-1} B_i^\top P_{i+1} A_i}_{\text{feedback gain}} \qquad \underbrace{\alpha_i = -(R_i + B_i^\top P_{i+1}B_i)^{-1} B_i^\top (P_{i+1}d_i + p_{i+1})}_{\text{affine (linearization-offset) term}} $$
$$ P_i = Q_i + K_i^\top R_i K_i + (A_i-B_iK_i)^\top P_{i+1}(A_i-B_iK_i), \qquad p_i = Q_ie_i - K_i^\top R_i\alpha_i + (A_i-B_iK_i)^\top\big(P_{i+1}(B_i\alpha_i+d_i) + p_{i+1}\big) $$

with $P_N$ the terminal Hessian. Backward pass costs $N$ solves of size $n_u$ (that inverse) plus a handful of $n_x$-sized matrix products; forward pass applies $\delta u_i = -K_i\,\delta x_i + \alpha_i$. Read the output as a controller, not a vector: the QP's minimizer is not just numbers $Z^\star$ but the affine feedback law $\kappa_i(\cdot) = \bar u_i - K_i(\cdot - \bar x_i) + \alpha_i$ — per-stage gains with bounded complexity. acados hands you exactly these ($K_i$) whenever you ask, which is the seed of §6.3. The interactive demo (§8) condenses these recursions into the box-QP and solves it with a tiny barrier IPM, so it is exactly this Riccati core that the QP wrapper iterates — the same layering HPIPM uses.

Numbers for scale. Cart-pole, $N=16$: full condensing gives a dense KKT of size $\approx (17\cdot4{+}16) = 84$ — a $84^3$ Cholesky is nothing. A 6-DOF robot arm, $n_x=12$, $N=40$: condensed KKT $\approx 600$, factorization $O(600^3)$ — expensive; Riccati is $40\times$ (a few $6\times6$ solves) — microseconds. The crossover is exactly the $n_x$-vs-$N$ balance in Fig. 5.
Takeaway: HPIPM gives acados three organizations of the same KKT system — Riccati (linear in $N$), partial condensing (tunable), dense (tiny problems) — chosen by problem geometry, all on BLASFEO panel-major linear algebra with no runtime allocation.

6The real-time iteration: warm-starting turns "solve an NLP" into "one QP"

Motivation. Sections 3–5 make each SQP iteration cheap; the RTI makes the count of iterations per sample equal to one. The enabling observation is control-specific: MPC problems are solved in a sequence, with a moving but slowly changing truth. The solution at $t$ is a near-perfect initial guess for $t{+}\Delta t$, because the new problem is the old one with one stage chopped off, one stage appended, and $\hat x_0$ replaced by a nearby measurement.

Figure 6 — the RTI cycle: shift, one Newton step, close the loop
plan at $t$: $x^{(t)}_{0},\dots,x^{(t)}_{N}$
warm-start guess at $t{+}\Delta t$: old plan shifted one stage
$\hat x_0^{(t)}$ applied, plant moves to $\hat x_0^{(t+\Delta t)}$
old tail re-used as guess; new stage appended
at $t{+}\Delta t$: ONE linearization @ shifted guess + ONE QP solve
apply $\Delta u_0$ + keep gains $K_i$ as feedback
total budget: "SQP with 1 iteration" — sub-ms with HPIPM+BLASFEO
shift guess
repeat every $\Delta t$, forever
The guess quality is excellent (one stage shorter horizon + small disturbance), so one Newton step lands near the new optimum. What you deliberately give up: a small residual optimality error $\|Z^\star - Z^{(t)}\|$, which theory bounds as proportional to the sampling period — see below.

6.1 Why "not converged" is not "wrong"

This looks like cheating, and the papers (Diehl, Lu, Schlöder — the origin of RTI; Frey/Diehl for acados' RTI-QP implementation) make it precise. Define the per-cycle loop as the map $\Phi_\Delta(Z^{(t)}, w_t)$ = "shift, one Newton step at sampling period $\Delta$". Classical Newton gives $\|Z^{(t+1)} - Z^\star_\Delta\| \le \rho\,\|Z^{(t)}-Z^\star_\Delta\| + c\,\Delta$ (contraction $\rho\lt 1$ near the solution, plus an $O(\Delta)$ term from the problem actually moving). Unrolling with zero warm-start error would give $O(\Delta/(1-\rho))$ — but with RTI, the steady-state distance to the true optimal solution is $O(\Delta)$: one iteration per cycle already achieves the same order of accuracy as full convergence, because the dominant error is the sampling disturbance, not the iteration. The iteration is linearly convergent with respect to the moving target. Iterating to tolerance with a stale measurement optimizes yesterday's problem more precisely; the extra accuracy buys nothing and costs the sample deadline.

6.2 The integrator's role in RTI: sensitivity propagation

One subtlety most explanations miss: at the first stage the state is now known ($\hat x_0$ replaced by a measurement), so the RTI linearization must be anchored there. acados' "SQP with 1 iteration + sensitivity propagation" does the Riccati pass and forwards the affine deviation rollout (the $d_i$-driven pass of §5.2) through the integrator sensitivities, which is exactly the GNSF/integrator design of §3 being asked to deliver cheap $\partial f_i$ per stage. That rollout is the trajectory the QP "predicts"; the applied input is the first element.

6.3 The feedback you get for free

Because the QP is solved by Riccati, its solution is a feedback law: the demo's applied input is $u_0 = \bar u_0 + \Delta u_0 + K_0(\hat x_0^{\rm meas} - \bar x_0)$ — the offset term of the affine policy evaluated at the actual measurement. Between samples, disturbances are damped by the Riccati gains (this is the "feedback part of the RTI", implemented in acados as the RTI update with the linearized model). A batch-solve dense-QP RTI can also return these gains (dual-solve / condensed-gain tricks), but Riccati produces them as a natural by-product — another vote for the block-tridiagonal path when $n_x$ is small.

Takeaway: RTI replaces "solve the NLP each sample" with "take one Newton step from a shifted warm start each sample" — provably $O(\Delta t)$-accurate for the moving problem, while the Riccati gains carried from the QP turn the plan into a feedback law. This is the exact schedule acados runs when you set nlp_solver = SQP_RTI (with sim_method_sens* options for sensitivity propagation, and the adjoint/GN gain update for the feedback part).

7The machine underneath: BLASFEO, codegen, and the solver stack

Motivation. Algorithmic structure sets the asymptotic work; embedded solvers live or die on constants. Two design decisions in acados' C core determine its wall-clock numbers, and both are worth understanding because they explain what acados is not (not a Python port of an NLP solver; not a fully code-generated monolith).

7.1 BLASFEO: why matrix–matrix products beat hand-written loops

Every Riccati step, condensing pass, and IPM linear solve is a short sequence of level-3 BLAS calls ($C \mathrel{+}= AB$, triangular solves). Ordinary BLAS is tuned for $1000^+$-dimension matrices; MPC wants dozens of $4{\times}4$ to $50{\times}50$ ones, where library call overhead, cache-unfriendly strides, and poor SIMD utilization dominate. BLASFEO (acados' linear-algebra layer, same authors) attacks exactly this regime with a panel-major matrix format: small dense blocks (panels) stored in a hardware-aware, SIMD-vector-friendly layout, plus an API that lets codegen emit fused calls for known dimensions. Published results: up to ~7× over optimized BLAS for the small matrices of embedded MPC, and it runs on x86, ARMv7/v8, PowerPC — the same binary compiles for a car ECU and a laptop.

7.2 What acados code-generates, and what it doesn't

Predecessor ACADO code-generated the entire solver (fixed size, fixed algorithm, self-contained C — fast but brittle). acados deliberately narrows codegen to the parts where it pays unambiguously:

Figure 7 — the acados stack: where the milliseconds go
Python / MATLAB–Simulink interface CasADi MX symbolic expressions for $\ell_i, f_c, h_i$; AcadosOcp description; JSON/C generation via acados_template
C codegen (the deliberate exception) per-function casadi_function-style evaluation of model/cost/constraints and their derivatives — generated straight-line C compiled into the solver's model module
Acados OCP solver (hand-written C, runtime-configurable) SQP loop · integrators (ERK/IRK/GNSF) · RTI shift & sensitivity propagation · line search · ocp_nlp_*_c
QP framework: HPIPM OCP-QP / dense-QP / tree-QP types · Riccati (RTPI) · partial/full condensing + expansion · primal-dual IPM variants · external qpOASES/OSQP/qpDUNES option
BLASFEO panel-major dense linear algebra, no runtime allocation, embedded-friendly

Why not codegen the solver? It was tried (ACADO) and traded flexibility for speed that BLASFEO now provides anyway: with level-3 BLAS in the inner loop, the surrounding control logic is not the bottleneck, so it stays ordinary maintainable C. You can swap the QP solver, integrator, or SQP variant at run time by setting an option, without regenerating anything. What must be generated is only the derivative evaluation — symbolic differentiation is a compile-time-only concern.

Memory discipline. acados/HPIPM take workspace pointers from the caller and never call malloc internally; problem dimensions are set once. That is the difference between a demo that runs and a controller with a certified worst-case cycle time.

Cost of one RTI cycle is then roughly: $N\times$(1 integrator stage set + 1 sensitivity set) + Riccati/condensing $O(Nn_x^3)$ + IPM iterations (a handful) × each $\approx$ Riccati — everything else is bookkeeping.

Read this as "where the time is": the green layer is the only auto-generated part; the algorithmic layers above BLASFEO are compiled once, reused for all problems.

7.3 The menu you actually configure

KnobOptionsPick it when…
nlp_solver_typeSQP_RTI, SQP, DDP, SQP_WITH_FEASIBLE_QP Online MPC → SQP_RTI (one step/cycle). Offline trajectory optimization → SQP iterated to tolerance. DDP variants for smooth input-free problems.
qp_solverPARTIAL_CONDENSING_HPIPM (default), FULL_CONDENSING_HPIPM, FULL_CONDENSING_QPOASES, PARTIAL_CONDENSING_OSQP, PARTIAL_CONDENSING_QPDUNES, … Small $n_x$, long $N$ → partial condensing with $N_c$ large, or the Riccati-IPM (RTI_RICATI_DP — what the C RTI layer selects automatically). Tiny problems → full condensing. Warm-started active set → qpOASES.
qp_solver_cond_N$1\ldots N$ (macro-stages $N_c$)The Riccati↔dense dial of Fig. 5. $N_c{=}N$ ⇒ effectively no condensing (pure Riccati); small $N_c$ ⇒ dense tail. Tune for density vs. structure.
sim_method / model.integrator_typeERK (order 1–8), IRK (Radau IIA), GNSF, LIFTED_IRK, DISCRETE; CVODES via C interface§3's decision tree.
tf / time_stepshorizon $T$, mesh refinement, $N$$N\Delta t$ fixed: bigger $\Delta t$ ⇒ fewer stages but more integrator work + accuracy loss per stage.
nlp_solver_tol_*, print_timeresidual tolerances, timingDiagnosing "not converging": residual $=$ max over $\|r_{\rm gap}\|$, dual, complementarity.
constraintsper-stage $h$, box, cost bound, soft + slack penalty, con_hessian_type Soft constraints (L1 or L2 slack) are the standard fix for infeasibility from disturbances.
MHE (acados_mhe_solver)arrival-cost recursion, same integratorsMoving-horizon estimation is the same machinery run on a different cost structure — see §9.
Takeaway: performance comes from three non-negotiable engineering choices — panel-major BLAS for small matrices (BLASFEO), codegen restricted to derivative evaluation, and allocation-free solver code — which together make the algorithmic structure of §3–§6 measurable rather than theoretical.

8Worked example, end to end: cart-pole balance and push recovery

Motivation. Everything above should now be traceable in one concrete cycle. Cart-pole, upright coordinates ($\theta{=}0$ is upright):

$$ \ddot p = \frac{u + m\sin\theta\,(l\dot\theta^2 - g\cos\theta)}{M + m\sin^2\theta}, \qquad \ddot\theta = \frac{g\sin\theta - \ddot p\cos\theta}{l}, $$

with $M{=}1$, $m{=}0.2$, $l{=}0.5$, $g{=}9.81$; $x=(p,\dot p,\theta,\dot\theta)$, $\theta$ measured from the upright; stage cost $(x-x^{\rm ref})^\top Q(x-x^{\rm ref}) + r u^2$, $Q=\mathrm{diag}(30,2,200,4)$, $r=0.02$; the terminal weight is the discrete DARE solution at the upright (the standard recipe that makes a finite horizon behave like an infinite-horizon LQR — computed by iterating the Riccati map of §5.2 to convergence, $P_N \approx \mathrm{diag}(628,236,1911,76)$ plus cross-terms); $|u|\le 25$ N enters the QP as a box constraint, not a post-hoc clip; $N=16$, $\Delta t = 0.05$ s. One RTI cycle at $t$:

  1. Shift the previous nominal $(\bar x_1,\bar u_1),\dots,(\bar x_N,\cdot)$ down one index, append an extrapolated final stage. Cost: a pointer/copy — $O(N n_x)$.
  2. Anchor $\bar x_0 \leftarrow \hat x_0^{\rm meas}$; $\delta x_0 = 0$ in deviation coordinates.
  3. Linearize all 16 stages: ERK-4 gives $(A_i,B_i)$ and the shooting-gap residual $r_i = \Phi(\bar x_i,\bar u_i) - \bar x_{i+1}$ (zero only if the nominal is feasible).
  4. Condense the linearized gaps into a dense QP over the $N{=}16$ input increments $\delta u$ with box constraints $|\bar u_i + \delta u_i|\le 25$ — the $O(N^2 n_x^3)$ build-up of §5, with the DARE-weighted cost (§5.2 recursion, wrapped angle error).
  5. Solve the box-QP with the tiny barrier interior-point method (a fixed handful of Newton iterations, warm-started from last cycle's solution — the HPIPM pattern).
  6. Roll forward $\delta x_{i+1} = A_i\delta x_i + B_i\delta u_i + r_i$; new nominal $x_i^+ = \bar x_i + \delta x_i$, $u_i^+ = \bar u_i + \delta u_i$; apply $u_0^+$. One cycle, one QP solve, done — at $n_x{=}4,n_u{=}1$ the box-QP is 16 scalar barrier-Newton iterations: microseconds in C.

The three regimes, live

The following widget runs the exact algorithm (Gauss–Newton multiple shooting + condensing + box-QP barrier solve + RTI shift/warm start) in JavaScript on a 30 Hz virtual clock, so you can feel the design decisions rather than read about them. Three settings, same plant, same weights. ⏸/▶ pauses and resumes the loop, “Push the pole” injects a $\dot\theta_0=5.5\ \mathrm{rad/s}$ kick, and the controller's measurement carries a fixed, seeded noise realization so every trace stays visibly alive (switch modes or press Reset to replay from the same start):

p (m)
θ (rad)
+1.0 … −1.0
u (N)
+25 … −25
shooting
gap
log scale
1e-10 … 1e3
t = 0.00 s · ready
Figure 8 — what to look for in the demo

RTI (warm start, 1 SQP step)

Balances and recovers from a hard push ($\dot\theta_0 = 5.5\ \mathrm{rad/s}$: the pole swings to $\sim0.8$ rad and the cart lurches to $\approx$2.7 m before the controller hauls both back to the reference in a couple of seconds). The gap panel shows the signature: at the push the shooting gap spikes to $\sim1$ (one Newton step can't fully close a big gap), then decays below $10^{-6}$ within a couple of seconds as the warm start re-converges — the loop never runs the 25-iteration inner loop, and it doesn't need to (§6.1's $O(\Delta t)$ bound).

Cold start, 1 SQP step

Identical work per cycle, no warm start: every cycle the linearization restarts from a flat nominal, so the single Newton step never accumulates trajectory knowledge. Mild operation works, but the same push makes the closed loop diverge — the direct evidence that RTI's speed is bought with warm starting, not with a magic solver.

Cold start, 25 SQP steps

Converges to the exact optimum every sample (gap at machine precision, $\sim10^{-16}$) — same closed-loop behavior as RTI, but $\approx$25× the SQP iterations ($\approx$100× the QP inner work) per sample. At 20 Hz on real hardware that is the difference between a controller and a space heater.
The log-gap panel is the honest scoreboard: RTI and full-convergence ride the same blue trajectory, while RTI spends ~25× fewer SQP steps (~100× less QP inner work) per cycle; its gap spikes only while a disturbance is being rejected. Cold-start-1 spends the same as RTI but has no plan to inherit, and pays for it the moment a disturbance arrives.
Takeaway: one RTI cycle = shift + measure + linearize + solve QP + apply; the gap panel makes the bargain visible — RTI's gap spikes while rejecting a disturbance, decays on its own, and the loop never needs the 25-iteration inner solve; only the warm start distinguishes it from a cold start that diverges under the same disturbance.

9Variations, siblings, and where this fits

9.1 Same machinery, other problems

Moving horizon estimation (MHE)

acados' MHE solver is the same SQP + integrator + QP stack on a different cost structure: arrival cost $V(\hat x_0)$ (propagated by a Riccati-like recursion, optionally the "covariance approximation") plus measured-output residuals over the window. Estimation and control then share one code path and one integrator — the reason acados is popular for NMPC+MHE pairs.

Tree-structured / scenario / robust MPC

HPIPM's third QP type replaces the stage chain with a tree ($x_m = A_nx_n + B_nu_n$ for children $m\in C(n)$). Scenario-based stochastic MPC and multistage robust MPC are the same Riccati/condensing idea with a branching index replacing the time index — a nice illustration of how general the "eliminate the coupling variables in a topological order" trick is.

Closed-loop NMPC: tube, LTV, learning

Because the QP hands back $K_i$, acados is the natural host for tube MPC (nominal LP/QP plan + tightened constraints + a pre-stabilizing $K$ absorbing the disturbance), LTV-MPC (linearize once, keep gains — the low-compute cousin of RTI), and learned-policy-with-safety-layer stacks where acados is the certificate.

Parametric NLP / MPC-RL hybrids

acados exposes sensitivities of the solution w.r.t. parameters ($\partial Z^\star/\partial p$), used for closed-loop sensitivity computation, MPC as a differentiable layer, and MPC-RL blending (policy proposes, acados projects onto feasible/stabilizing inputs).

9.2 Where acados sits among its cousins

FamilyExamplesStructural difference from acados
Code-generated SQP/RTIACADO, nmpc-codegen, FORCESpro, AutoGenU Whole solver specialized per problem: fastest constants, but fixed dimensions/algorithm and large build times. acados keeps the solver generic on top of BLASFEO, and only codegens derivatives.
Generic sparse-NLP / IPMIPOPT, CasADi+IPOPT, KNITRO, SOCP-based No multi-stage assumption: a generic sparse KKT, superlinear convergence per iteration (needs several iterations), no Riccati, no RTI. Right for offline trajectory optimization; wrong for a 1 ms loop.
Convex-only MPC (LTV/LMPC)OSQP / qpOASES alone, HPMPC, muAO-MPC, DuQuad Solve one QP per sample, no nonlinear model. Fast and provable — but only for linear (or once-linearized) plants.
Shooting-free / collocationDIRCOL, GPOPS Full discretization at many collocation points, huge sparse NLP: great for one high-accuracy trajectory, not for repeated online solutions.

9.3 The honest limits

The design has prices worth naming. (i) One QP per sample means the warm start must be good: large measurement jumps, mode switches, or a contact event can leave the RTI step far from the new optimum — mitigations are more SQP iterations (feasibility_only steps), soft constraints, or re-linearizing from the measured state (SQP with a few steps). (ii) First-order model accuracy is your accuracy: GN ignores $\partial^2 f$ curvature, so a coarse mesh or a badly conditioned cost model shows up as steady-state offset or slow convergence. (iii) The integrator is on the critical path: stiff plants force IRK/CVODES, whose per-stage cost can dominate the Riccati work, and CVODES' variable work breaks worst-case timing. (iv) Stability is not free: RTI is $O(\Delta t)$-accurate but a stabilizing NMPC needs the standard ingredients — terminal cost/region (Lyapunov conditions), sufficiently small $\Delta t$, or contraction-based conditions; the solver's speed does not create these. acados gives you the milliseconds to enforce them inside the horizon, nothing more.

Takeaway: acados = multiple shooting + Gauss–Newton SQP + structure-exploiting QP (Riccati/condensing in HPIPM) + real-time iteration, all on BLASFEO. Every layer exists to convert a specific piece of control-problem structure into wall-clock time — and the whole stack is roughly "Newton on the KKT system of a multiple-shooting discretization, done with the structure left intact."

10References and further reading