B-Spline-Based Model Predictive Control

• B-spline-based MPC is a receding-horizon control framework that parameterizes continuous trajectories via spline coefficients, ensuring smooth derivatives.
• It leverages differential flatness to decouple decision variables from collocation points, achieving significant computational speedups over dense direct transcription methods.
• Variants using Bézier curves and cubic Hermite splines extend the approach to convex, sampling-based, and adaptive predictive control, enabling robust constraint handling.

Searching arXiv for recent and relevant papers on B-spline-based MPC to ground the article and verify related work. Searching arXiv for "B-spline model predictive control", "Bézier MPC", and "flatness-based MPC B-splines". B-spline-based Model Predictive Control (MPC) denotes a class of receding-horizon control methods in which the predicted state, input, or flat-output trajectories are parameterized by spline basis functions rather than by a dense sequence of time-sampled decision variables. In the pusher-slider setting studied by Neve et al., the combination of differential flatness with a B-splines transcription replaces a large direct-transcription nonlinear program by a compact optimization over spline control points, yielding empirical computational acceleration up to $65\%$ relative to a Direct Multiple Shooting baseline (Neve et al., 2023). Across adjacent strands of the literature, closely related spline parameterizations include Bézier curves for convex multi-rate MPC and control Lyapunov function tracking (Csomay-Shanklin et al., 2022), dual-space cubic Hermite spline parameterizations within sampling-based MPC (Schramm et al., 24 Nov 2025), and cubic-Hermite B-spline observables in adaptive predictive control (Alhazmi et al., 5 Feb 2026). Taken together, these works situate B-spline-based MPC as a family of trajectory-parameterized predictive controllers rather than a single algorithmic template.

1. Conceptual basis and relation to direct transcription

In standard MPC, at each sampling instant one solves a finite-horizon optimal control problem

$$\text{minimize over }u(\cdot):\quad J=\phi(x(T))+\int_0^T \ell(x(t),u(t))\,dt$$

subject to

$$\dot x(t)=f(x(t),u(t)),\qquad x(0)=x_0,\qquad x(T)=x_{\text{goal}},\qquad h(x(t),u(t))\le 0.$$

After solving, the first portion of the optimal control $u^*(t)$ is applied, and at the next sampling instant the horizon is shifted, the plant is re-measured, and a new optimal control problem is solved (Neve et al., 2023).

The classical alternative emphasized in the pusher-slider work is “direct” transcription, exemplified by Direct Multiple Shooting. There, the horizon is discretized into $N$ intervals, one introduces both states $x_0,\dots,x_N$ and controls $u_0,\dots,u_{N-1}$ as optimization variables, and the dynamics are enforced through equality constraints

$x_{k+1}-\Phi(x_k,u_k)=0,\qquad k=0,\dots,N-1,$

where $\Phi$ is a discrete-time integration map such as Runge–Kutta. The resulting nonlinear program typically has $O(N\cdot(n_x+n_u))$ variables and $$\text{minimize over }u(\cdot):\quad J=\phi(x(T))+\int_0^T \ell(x(t),u(t))\,dt$$0 equality constraints, so for real-time applications with sampling times in the tens of milliseconds, online solution can be prohibitive (Neve et al., 2023).

B-spline-based MPC changes the discretization philosophy. Instead of optimizing every state-control sample directly, it parameterizes a continuous trajectory by a comparatively small set of spline coefficients. In the flatness-based pusher-slider formulation, the decision variable is the vector of spline control points $$\text{minimize over }u(\cdot):\quad J=\phi(x(T))+\int_0^T \ell(x(t),u(t))\,dt$$1, with $$\text{minimize over }u(\cdot):\quad J=\phi(x(T))+\int_0^T \ell(x(t),u(t))\,dt$$2 because each $$\text{minimize over }u(\cdot):\quad J=\phi(x(T))+\int_0^T \ell(x(t),u(t))\,dt$$3 (Neve et al., 2023). This suggests that the principal computational benefit comes not from eliminating constraints altogether, but from decoupling the number of optimization variables from the number of collocation points.

2. B-spline trajectory parameterization and flatness-based lifting

For the quasi-static pusher-slider, one may choose the Cartesian coordinates of the slider’s center $$\text{minimize over }u(\cdot):\quad J=\phi(x(T))+\int_0^T \ell(x(t),u(t))\,dt$$4 as flat outputs. Differential flatness means that all states $$\text{minimize over }u(\cdot):\quad J=\phi(x(T))+\int_0^T \ell(x(t),u(t))\,dt$$5 and inputs $$\text{minimize over }u(\cdot):\quad J=\phi(x(T))+\int_0^T \ell(x(t),u(t))\,dt$$6 can be written as algebraic functions of a flat output $$\text{minimize over }u(\cdot):\quad J=\phi(x(T))+\int_0^T \ell(x(t),u(t))\,dt$$7 and a finite number of its derivatives: $$\text{minimize over }u(\cdot):\quad J=\phi(x(T))+\int_0^T \ell(x(t),u(t))\,dt$$8 For the pusher-slider, the remaining state $$\text{minimize over }u(\cdot):\quad J=\phi(x(T))+\int_0^T \ell(x(t),u(t))\,dt$$9 and orientation $$\dot x(t)=f(x(t),u(t)),\qquad x(0)=x_0,\qquad x(T)=x_{\text{goal}},\qquad h(x(t),u(t))\le 0.$$0, as well as the tangential and normal push velocities $$\dot x(t)=f(x(t),u(t)),\qquad x(0)=x_0,\qquad x(T)=x_{\text{goal}},\qquad h(x(t),u(t))\le 0.$$1, can be expressed via derivatives of $$\dot x(t)=f(x(t),u(t)),\qquad x(0)=x_0,\qquad x(T)=x_{\text{goal}},\qquad h(x(t),u(t))\le 0.$$2 and $$\dot x(t)=f(x(t),u(t)),\qquad x(0)=x_0,\qquad x(T)=x_{\text{goal}},\qquad h(x(t),u(t))\le 0.$$3 (Neve et al., 2023): $$\dot x(t)=f(x(t),u(t)),\qquad x(0)=x_0,\qquad x(T)=x_{\text{goal}},\qquad h(x(t),u(t))\le 0.$$4

$$\dot x(t)=f(x(t),u(t)),\qquad x(0)=x_0,\qquad x(T)=x_{\text{goal}},\qquad h(x(t),u(t))\le 0.$$5

$$\dot x(t)=f(x(t),u(t)),\qquad x(0)=x_0,\qquad x(T)=x_{\text{goal}},\qquad h(x(t),u(t))\le 0.$$6

while $$\dot x(t)=f(x(t),u(t)),\qquad x(0)=x_0,\qquad x(T)=x_{\text{goal}},\qquad h(x(t),u(t))\le 0.$$7 depends on third-order derivatives. Because of flatness, any sufficiently smooth $$\dot x(t)=f(x(t),u(t)),\qquad x(0)=x_0,\qquad x(T)=x_{\text{goal}},\qquad h(x(t),u(t))\le 0.$$8 uniquely defines $$\dot x(t)=f(x(t),u(t)),\qquad x(0)=x_0,\qquad x(T)=x_{\text{goal}},\qquad h(x(t),u(t))\le 0.$$9 without additional inter-node continuity constraints (Neve et al., 2023).

The spline transcription begins by choosing a nondecreasing knot vector

$u^*(t)$0

where $u^*(t)$1 is the spline degree. For a clamped uniform B-spline, one sets $u^*(t)$2, $u^*(t)$3, and places the interior knots equispaced in $u^*(t)$4. The B-spline basis functions $u^*(t)$5 are defined by the Cox–de Boor recursion (Neve et al., 2023): $u^*(t)$6

$u^*(t)$7

The flat trajectory is then parameterized over $u^*(t)$8 by $u^*(t)$9 control points $N$0: $N$1 Hence

$N$2

The derivatives $N$3 are obtained by differentiating the B-splines, which remain B-splines of degree $N$4 (Neve et al., 2023). At any time $N$5, once $N$6 are known, one computes

$N$7

$N$8

The entire state-input trajectory is therefore a deterministic function of the control-point vector $N$9 (Neve et al., 2023).

A closely related construction appears in the multi-rate framework of “Multi-Rate Planning and Control of Uncertain Nonlinear Systems: Model Predictive Control and Control Lyapunov Functions,” which uses Bézier curves, explicitly described there as a special case of B-spline, to stitch together $x_0,\dots,x_N$0 segments of order $x_0,\dots,x_N$1 over subintervals of length $x_0,\dots,x_N$2 (Csomay-Shanklin et al., 2022). In that setting, adjoining segments match in value and derivative up to order $x_0,\dots,x_N$3, and each Bézier segment lies in the convex hull of its control points. The common structural point is that spline control points parameterize a continuous admissible trajectory while moving continuity and smoothness properties into the basis itself (Csomay-Shanklin et al., 2022).

3. Finite-dimensional optimization forms

Once the trajectory is parameterized by spline coefficients, the infinite-dimensional optimal control problem becomes a finite-dimensional optimization problem defined over those coefficients. In the pusher-slider transcription, one chooses collocation times $x_0,\dots,x_N$4, such as Greville points or uniformly spaced points in $x_0,\dots,x_N$5, and enforces boundary and path constraints there. Boundary conditions include

$x_0,\dots,x_N$6

and, if needed, $x_0,\dots,x_N$7 and $x_0,\dots,x_N$8. Path constraints are imposed on state and input: $x_0,\dots,x_N$9 Since the flatness relations eliminate the need for dynamic-continuity constraints, the only equality constraints are the boundary conditions plus any path equality constraints. When written in affine form in the control-point vector $u_0,\dots,u_{N-1}$0, the constraints become

$u_0,\dots,u_{N-1}$1

Here, $u_0,\dots,u_{N-1}$2 and $u_0,\dots,u_{N-1}$3 collect combinations of $u_0,\dots,u_{N-1}$4 and their derivatives (Neve et al., 2023).

The cost can be expressed directly in the spline coefficients. In the flatness-based pusher-slider case,

$u_0,\dots,u_{N-1}$5

Each $u_0,\dots,u_{N-1}$6 and $u_0,\dots,u_{N-1}$7 is a nonlinear but smooth function of $u_0,\dots,u_{N-1}$8. One can form a Gauss–Newton approximation or use a full-SQP solver; if $u_0,\dots,u_{N-1}$9 and $x_{k+1}-\Phi(x_k,u_k)=0,\qquad k=0,\dots,N-1,$0 are quadratic and the flatness maps are linearized, one obtains a local quadratic “QP-like” subproblem in $x_{k+1}-\Phi(x_k,u_k)=0,\qquad k=0,\dots,N-1,$1 (Neve et al., 2023). The resulting compact formulation is

$x_{k+1}-\Phi(x_k,u_k)=0,\qquad k=0,\dots,N-1,$2

The multi-rate Bézier framework yields a different optimization structure. There, the planner operates over discrete waypoints $x_{k+1}-\Phi(x_k,u_k)=0,\qquad k=0,\dots,N-1,$3 with $x_{k+1}-\Phi(x_k,u_k)=0,\qquad k=0,\dots,N-1,$4, while the Bézier control points are implicitly determined by adjacent waypoint pairs via the linear relation

$x_{k+1}-\Phi(x_k,u_k)=0,\qquad k=0,\dots,N-1,$5

State constraints are handled by convex-hull tightening over the control points, and input constraints are cast as second-order cone constraints involving auxiliary slacks $x_{k+1}-\Phi(x_k,u_k)=0,\qquad k=0,\dots,N-1,$6. If the stage and terminal costs are quadratic, the entire problem is an SOCP (Csomay-Shanklin et al., 2022).

The literature therefore contains multiple finite-dimensional optimization archetypes under the broader heading of spline-based MPC: compact nonlinear programs over B-spline coefficients (Neve et al., 2023), convex SOCPs over waypoints with Bézier lifting (Csomay-Shanklin et al., 2022), and, in a different computational regime, sampling-based path-integral optimization over cubic Hermite control points (Schramm et al., 24 Nov 2025). A plausible implication is that the spline representation itself is orthogonal to the final solver class.

4. Constraint handling, smoothness, and continuity

A central technical advantage of spline parameterizations is that continuity is encoded by basis construction rather than by explicit continuity constraints between all discretization nodes. In the pusher-slider formulation, flatness plus B-spline parameterization means that any sufficiently smooth flat output $x_{k+1}-\Phi(x_k,u_k)=0,\qquad k=0,\dots,N-1,$7 induces a unique state-input trajectory, and “the need for dynamic-continuity constraints” is eliminated (Neve et al., 2023). This shifts the constrained optimization burden toward endpoint conditions and collocated path constraints.

Smoothness requirements can be stringent. In the pusher-slider system, $x_{k+1}-\Phi(x_k,u_k)=0,\qquad k=0,\dots,N-1,$8 depends on third-order derivatives, so the implementation chooses B-spline degree $x_{k+1}-\Phi(x_k,u_k)=0,\qquad k=0,\dots,N-1,$9 “to ensure continuity up to 4th derivatives, needed for $\Phi$0” (Neve et al., 2023). This is an instance of a more general design principle: spline degree is determined not only by interpolation quality but by the derivative order required in state or input reconstruction.

In the multi-rate Bézier framework, constraint handling exploits the convex-hull property. Each Bézier segment $\Phi$1 lies in the convex hull of its $\Phi$2 control points $\Phi$3, and if those control points lie in the tightened polytope $\Phi$4, then the nonlinear state remains in $\Phi$5 for all time under the low-level CLF tracker (Csomay-Shanklin et al., 2022). Input limits are enforced through second-order cone constraints derived from bounds on the low-level feedback law. The paper states that the planned tube $\Phi$6 remains in $\Phi$7, the CLF tracker never violates $\Phi$8, and the shifted previous solution remains feasible at every replan under the theorem’s assumptions (Csomay-Shanklin et al., 2022).

The sampling-based work on reference-free locomotion uses a dual-space spline parameterization with position and velocity control points $\Phi$9, reconstructing

$O(N\cdot(n_x+n_u))$0

$O(N\cdot(n_x+n_u))$1

There, “intrinsic smoothness” is highlighted explicitly: by construction, any Hermite-interpolated trajectory is continuously differentiable, which cuts down wasted samples that would otherwise be spiky or physically infeasible (Schramm et al., 24 Nov 2025). Although that framework is not posed as a classical deterministic NLP, it reinforces the same underlying point: spline bases regularize admissible trajectory shape before optimization begins.

5. Computational scaling and real-time operation

The best-documented computational comparison in the supplied literature is the pusher-slider case study. Direct Multiple Shooting with $O(N\cdot(n_x+n_u))$2 shooting nodes uses $O(N\cdot(n_x+n_u))$3 decision variables plus $O(N\cdot(n_x+n_u))$4 equality constraints, whereas B-spline MPC uses only $O(N\cdot(n_x+n_u))$5 variables, independent of the number of collocation points $O(N\cdot(n_x+n_u))$6; only the number of constraints grows with $O(N\cdot(n_x+n_u))$7 (Neve et al., 2023). Typical experiments choose $O(N\cdot(n_x+n_u))$8, such as $O(N\cdot(n_x+n_u))$9 and $$\text{minimize over }u(\cdot):\quad J=\phi(x(T))+\int_0^T \ell(x(t),u(t))\,dt$$00 (Neve et al., 2023). In simulations on the pusher-slider, Neve et al. report “up to 65% reduction in average solve time per cycle” compared to a standard DMS-based MPC using IPOPT, attributing the speedup to the reduced variable count and sparser Jacobians and Hessians (Neve et al., 2023).

The reported real-time configuration in that work uses a horizon $$\text{minimize over }u(\cdot):\quad J=\phi(x(T))+\int_0^T \ell(x(t),u(t))\,dt$$01, $$\text{minimize over }u(\cdot):\quad J=\phi(x(T))+\int_0^T \ell(x(t),u(t))\,dt$$02 collocation points with $$\text{minimize over }u(\cdot):\quad J=\phi(x(T))+\int_0^T \ell(x(t),u(t))\,dt$$03, spline degree $$\text{minimize over }u(\cdot):\quad J=\phi(x(T))+\int_0^T \ell(x(t),u(t))\,dt$$04, and $$\text{minimize over }u(\cdot):\quad J=\phi(x(T))+\int_0^T \ell(x(t),u(t))\,dt$$05 control points with a clamped uniform knot vector with $$\text{minimize over }u(\cdot):\quad J=\phi(x(T))+\int_0^T \ell(x(t),u(t))\,dt$$06 knots (Neve et al., 2023). The solver is IPOPT called via CasADi, warm-starting from the previous $$\text{minimize over }u(\cdot):\quad J=\phi(x(T))+\int_0^T \ell(x(t),u(t))\,dt$$07 for faster convergence. A shrinking-horizon strategy is used: by re-parameterizing the continuous spline on $$\text{minimize over }u(\cdot):\quad J=\phi(x(T))+\int_0^T \ell(x(t),u(t))\,dt$$08 each cycle, one avoids rebuilding continuity constraints; only $$\text{minimize over }u(\cdot):\quad J=\phi(x(T))+\int_0^T \ell(x(t),u(t))\,dt$$09 change at the endpoints (Neve et al., 2023).

The convex Bézier-based multi-rate architecture emphasizes a different computational profile. Its constraints are all convex—linear dynamics, linear inequalities on the Bézier control points, and second-order cone constraints for input—and the paper notes that one obtains an SOCP of modest size, solvable by MOSEK, ECOS, and similar solvers; the simulations use MOSEK via Yalmip (Csomay-Shanklin et al., 2022).

By contrast, the reference-free sampling-based work argues for sample efficiency and CPU real-time operation through low-dimensional spline parameterization. Instead of sampling a high-frequency sequence of torques, it samples only $$\text{minimize over }u(\cdot):\quad J=\phi(x(T))+\int_0^T \ell(x(t),u(t))\,dt$$10 nodes per degree of freedom, reducing the number of random variables from $$\text{minimize over }u(\cdot):\quad J=\phi(x(T))+\int_0^T \ell(x(t),u(t))\,dt$$11 down to $$\text{minimize over }u(\cdot):\quad J=\phi(x(T))+\int_0^T \ell(x(t),u(t))\,dt$$12. The paper states that with only $$\text{minimize over }u(\cdot):\quad J=\phi(x(T))+\int_0^T \ell(x(t),u(t))\,dt$$13–$$\text{minimize over }u(\cdot):\quad J=\phi(x(T))+\int_0^T \ell(x(t),u(t))\,dt$$14 rollouts per $$\text{minimize over }u(\cdot):\quad J=\phi(x(T))+\int_0^T \ell(x(t),u(t))\,dt$$15–$$\text{minimize over }u(\cdot):\quad J=\phi(x(T))+\int_0^T \ell(x(t),u(t))\,dt$$16 ms on a multicore CPU, it achieves $$\text{minimize over }u(\cdot):\quad J=\phi(x(T))+\int_0^T \ell(x(t),u(t))\,dt$$17 Hz MPC without a GPU (Schramm et al., 24 Nov 2025). The computational message common to all three formulations is not a single complexity theorem, but a recurring dimensionality argument: spline coefficients can be substantially fewer than dense time-grid decision variables.

6. Variants, neighboring formulations, and interpretive issues

B-spline-based MPC is not restricted to flatness-based nonlinear programming. One neighboring formulation uses Bézier curves, described explicitly as “a special case of B-spline,” to enable planning continuous trajectories respecting constraints by planning a sequence of discrete points (Csomay-Shanklin et al., 2022). Another uses cubic Hermite interpolation; the locomotion paper states that a cubic Hermite interpolant can also be written “in the standard B-spline language by introducing an appropriate knot vector and repeated knots at the boundaries,” though the dual-space design with explicit velocity control points is easier to impose directly in the Hermite basis (Schramm et al., 24 Nov 2025). These works indicate that the term “B-spline-based MPC” often encompasses a broader spline-parameterized design space rather than only uniform B-splines in the Cox–de Boor form.

A related but distinct use of splines appears in adaptive predictive control. “Nonlinear Predictive Cost Adaptive Control of Pseudo-Linear Input-Output Models Using Polynomial, Fourier, and Cubic Spline Observables” employs cubic-Hermite B-spline observables in both online identification and receding-horizon optimization (Alhazmi et al., 5 Feb 2026). There, the scalar argument $$\text{minimize over }u(\cdot):\quad J=\phi(x(T))+\int_0^T \ell(x(t),u(t))\,dt$$18 is divided into $$\text{minimize over }u(\cdot):\quad J=\phi(x(T))+\int_0^T \ell(x(t),u(t))\,dt$$19 equal subintervals, and at each internal knot $$\text{minimize over }u(\cdot):\quad J=\phi(x(T))+\int_0^T \ell(x(t),u(t))\,dt$$20 two local cubic-Hermite shape functions $$\text{minimize over }u(\cdot):\quad J=\phi(x(T))+\int_0^T \ell(x(t),u(t))\,dt$$21 and $$\text{minimize over }u(\cdot):\quad J=\phi(x(T))+\int_0^T \ell(x(t),u(t))\,dt$$22 define the two-component basis

$$\text{minimize over }u(\cdot):\quad J=\phi(x(T))+\int_0^T \ell(x(t),u(t))\,dt$$23

which is stacked into the $$\text{minimize over }u(\cdot):\quad J=\phi(x(T))+\int_0^T \ell(x(t),u(t))\,dt$$24-dimensional spline basis

$$\text{minimize over }u(\cdot):\quad J=\phi(x(T))+\int_0^T \ell(x(t),u(t))\,dt$$25

This basis becomes the observable $$\text{minimize over }u(\cdot):\quad J=\phi(x(T))+\int_0^T \ell(x(t),u(t))\,dt$$26 in a pseudo-linear input-output model $$\text{minimize over }u(\cdot):\quad J=\phi(x(T))+\int_0^T \ell(x(t),u(t))\,dt$$27, with parameters updated by SIFt-RLS and then used within iterative MPC (Alhazmi et al., 5 Feb 2026). In Example 4, the cubic Hermite spline controller with four segments, denoted CB$$\text{minimize over }u(\cdot):\quad J=\phi(x(T))+\int_0^T \ell(x(t),u(t))\,dt$$28, achieved the smallest command-following error $$\text{minimize over }u(\cdot):\quad J=\phi(x(T))+\int_0^T \ell(x(t),u(t))\,dt$$29 and the lowest one-step prediction error across $$\text{minimize over }u(\cdot):\quad J=\phi(x(T))+\int_0^T \ell(x(t),u(t))\,dt$$30, and in Examples 5 and 6 it again outperformed the polynomial and Fourier alternatives listed in the paper (Alhazmi et al., 5 Feb 2026).

Several misconceptions are therefore best avoided. First, spline-based MPC is not inherently convex: the pusher-slider formulation is a compact nonlinear program (Neve et al., 2023), the multi-rate Bézier formulation is an SOCP under quadratic costs (Csomay-Shanklin et al., 2022), and the reference-free locomotion approach is zero-order MPPI rather than gradient-based optimization (Schramm et al., 24 Nov 2025). Second, spline use does not imply that only outputs are parameterized; depending on the formulation, the coefficients may represent flat outputs (Neve et al., 2023), stitched state trajectories (Csomay-Shanklin et al., 2022), joint positions and velocities simultaneously (Schramm et al., 24 Nov 2025), or nonlinear observables for system identification (Alhazmi et al., 5 Feb 2026). Third, computational savings arise from reduced decision dimension and structured smoothness, not from any claim that constraints disappear; in most formulations, constraints remain central and are transferred to control points, collocation points, tightened sets, or sampled trajectories.

Taken together, the cited works portray B-spline-based MPC as a methodological pattern with three recurring ingredients: low-dimensional continuous trajectory parameterization, derivative-structured smoothness inherited from the basis, and a receding-horizon optimizer whose burden shifts from dense node-by-node dynamics enforcement to a smaller set of spline coefficients. In the pusher-slider example this pattern enabled real-time MPC with documented $$\text{minimize over }u(\cdot):\quad J=\phi(x(T))+\int_0^T \ell(x(t),u(t))\,dt$$31 empirical speedup over a Direct Multiple Shooting baseline (Neve et al., 2023). In adjacent formulations it supports convex constraint handling through the convex-hull property (Csomay-Shanklin et al., 2022), sampling-based real-time locomotion without predefined contact schedules (Schramm et al., 24 Nov 2025), and adaptive pseudo-linear predictive control through cubic-spline observables (Alhazmi et al., 5 Feb 2026).