Exploiting Differential Flatness for Efficient Learning-based Model Predictive Control of Constrained Multi-Input Control Affine Systems
arXiv:2604.24706v1 [eess.SY] 27 Apr 2026
Tobias A. Farger1,∗ , Adam W. Hall2,∗ , and Angela P. Schoellig1 Abstract— Learning-based control techniques use data from past trajectories to control systems with uncertain dynamics. However, learning-based controllers are often computationally inefficient, limiting their practicality. To address this limitation, we propose a learning-based controller that exploits differential flatness, a property of many robotic systems. Recent research on using flatness for learning-based control either is limited in that it (i) ignores input constraints, (ii) applies only to single-input systems, or (iii) is tailored to specific platforms. In contrast, our approach uses a system extension and block-diagonal cost formulation to control general multi-input, nonlinear, affine systems. Furthermore, it satisfies input and half-space flat state constraints and guarantees probabilistic Lyapunov decrease using only two sequential convex optimizations. We show that our approach performs similarly to, but is multiple times more efficient than, a Gaussian process model predictive controller in simulation, and achieves competitive tracking in real hardware experiments.
I. I NTRODUCTION Learning-based control methods use machine learning models to estimate uncertain dynamics using data collected from previous trajectories. Recently, such methods have proven to be a less conservative, but more computationally expensive, option for controlling uncertain systems as compared to classical robust control [1], particularly when accurate first-principles models are difficult to obtain or system parameters vary across operating conditions. A common approach uses Gaussian processes (GPs) to model the dynamics within a model predictive control (MPC) framework [2], which necessitates solving a nonlinear program at every time step—making real-time computation challenging. One can, however, use structural assumptions about the system dynamics, like differential flatness [3], to improve the efficiency of learning-based controllers. Differential flatness, a property of many robotic systems (e.g., quadrotors [4] and flexible-joint manipulators [5]), allows for linear control techniques to be applied to nonlinear systems. This is because differentially flat systems can be exactly transformed into linear systems via nonlinear input transformations. Thus, linear control methods, like linear MPC, can be used to control flat nonlinear systems. Such flat model predictive control (FMPC) methods perform similarly to nonlinear model predictive control (NMPC) but require less compute [4], [6]. FMPC, however, cannot enforce input constraints while maintaining convexity and is sensitive to model mismatch in the nonlinear transformation. ∗ These authors contributed equally to the manuscript. 1 Learning Systems and Robotics Lab at the Technical University of Munich, Munich, Germany. 2 Learning Systems and Robotics Lab at the University of Toronto Institute for Aerospace Studies (UTIAS) and the Vector Institute for Artificial Intelligence, Toronto, Canada Email: {tobias.farger, angela.schoellig}@tum.de
[email protected] This work was supported by the Robotics Institute Germany under BMBF grant 16ME0997K.
Learned Flat Input Transformation v = Ψ(z, ū) GP 1 . . . GP m zkref
Flat MPC
vk∗ zk∗
µ, σ 2 Safety Filter as SOCP
ūk
System Extension
uk ẑk−1
Uncertain Differentially Flat System
1 Fig. 1. Our approach consists of a flat model predictive controller (MPC), a safety filter and a system extension. The flat MPC tracks the reference zkref and computes the optimal flat input vk∗ and state zk∗ . As the nonlinear flat input transformation v = Ψ(z, ū) is uncertain, it is modeled with an independent Gaussian processes (GPs) for each component. The second-order cone program (SOCP) safety filter then uses the GPs’ mean and covariance predictions µ and σ 2 , along with zk∗ and vk∗ , to find the extended input ūk that guarantees probabilistic Lyapunov decrease and flat state constraint satisfaction while respecting constraints on the input uk .
To eliminate any model mismatch and incorporate input constraints, recent work has explored learning-based controllers that exploit differential flatness. In [7], a robust linear quadratic regulator is combined with a learned nonlinear input mapping and its inverse for single-input systems. While this approach provides an upper bound on tracking error, it lacks the ability to explicitly handle input or flat state constraints. In contrast, [8] and [9] introduce a method tailored to multi-rotor systems, capable of handling multiple inputs and respecting input constraints. Although these techniques are efficient, they are limited to multi-rotor systems. Other recent approaches combine FMPC with a second order cone program (SOCP) safety filter that uses GPs to learn the nonlinear input transformation [10], [11]. This safety filter complements the FMPC by guaranteeing probabilistic Lyapunov decrease, probabilistic half-space flat state constraint satisfaction, and bounded input constraint satisfaction. However, these approaches are limited to single-input systems, and their input transformations do not require input derivatives (i.e., for systems with full relative degree [5]). The contributions of this work are: • A reformulated SOCP safety filter that supports general multi-input, control-affine systems with flat input transformations requiring input derivatives; • A novel Lipschitz-based bound and block-diagonal cost structure enabling probabilistic Lyapunov decrease; • Simulation and hardware validation on quadrotor trajectory tracking, showing competitive performance with GPMPC and NMPC, but with significantly reduced computation time. In contrast to prior work limited to single-input settings or quadrotor-specific formulations, we demonstrate a learningbased, flatness-exploiting MPC that scales to general multiinput systems while satisfying flat-state and input constraints with real-time computation. As in [11], our formulation still ensures probabilistic Lyapunov decrease and satisfaction of half-space flat state and input constraints for the next time step. Although the formulation does not guarantee recursive
feasibility, we show through simulation and hardware experiments that the controller remains feasible, stable, and constraint-satisfying in practice. Providing such guarantees for computationally efficient learning-based controllers is critical for deploying robotic systems that must operate in real time under uncertain dynamics and safety constraints. II. P ROBLEM S TATEMENT Consider a nonlinear multi-input control-affine system, ẋ(t) = f (x(t)) + g(x(t))u, (1) with initial condition x(0) = x0 . Here x(t) ∈ Rnx and u(t) ∈ Rm denote the system state and input, respectively. The unknown functions f : Rnx → Rnx and g : Rnx → Rnx ×m are assumed to be sufficiently smooth and (1) is known to be differentially flat with a known flat output. Definition II.1 (Differential Flatness [12]). A system (1) is differentially flat with respect to a flat output y = Λ(x, u, u̇, ..., u(ℓ) ) if all system states and inputs can be written as x = Φ(y, ẏ, ..., y (ρx ) ) and u = Γ(y, ẏ, ..., y (ρx +1) ). The input u is subject to constraints umin ≤ u ≤ umax (all vector inequalities are evaluated component-wise). In addition, we consider convex half-space constraints on the flat state z ∈ Z = {z|hT j z ≤ bj , ∀j = 1 . . . nc } ⊂ (ρ1 −1) nz , ..., m y, ..., m y (ρm −1) ]T , R , where P z = [1 y, ..., 1 y and nz = i (ρi − 1). We also have that ρx = maxi (ρi − 1). The subscript i at the bottom left indicates component i of a vector and (·)(p) denotes the derivative of order p. Remark II.1. Although we only consider constraints on the flat states z, these are often useful quantities like position, velocity, acceleration, and attitude. The objective is to design a computationally efficient controller for (1) that ensures accurate tracking of a flat state reference z ref and satisfies input and flat state constraints. III. BACKGROUND A. Discretized Exact Linear Dynamics Differential flatness implies exactly linear dynamics in a transformed flat-state space, allowing linear controllers to be designed in the flat space and mapped back to the original nonlinear system. Lemma III.1 (Linearized Flat Dynamics [5], [12]). Every flat system can be represented as a linear system and a nonlinear transformation ż = Az + Bv (2) v = Ψ(z, u, u̇, ..., u(ℓ) ), (ρ1 )
(ρm ) T
m
(3)
where v = [1 y , . . . , my ] ∈ R is the flat input, and A and B are such that (2) is a system of chained integrators. In this work, we discretize the flat system (2) such that zk+1 = Ad zk + Bd vk , (4) where k is the discrete time index, zk is the discretized flat state, vk is the discretized flat input, and Ad and Bd are exact discretizations of A and B, respectively. Since (2) decomposes into m independent integrator chains, Ad is block diagonal and Bd = diag(bd,1 , ..., bd,m ) with bd,i ∈ Rρi . Although this discretization introduces a sampling error, we assume that, if the sampling time is sufficiently small, this error can be easily accounted for by constraint tightening. Such discretizations are standard in the flat control literature [4], [6]–[9], [11].
B. Gaussian Processes We use GPs to model distributions over unknown but observable functions Ψ(a) : Rdim(a) → R. Prior knowledge about Ψ(a) is encoded via a prior mean and parameterized kernel function. Then, given a noisy data set Ψ̂(a) = Ψ(a)+ ε with ε ∼ N (0, σn2 ), the GP’s hyperparameters are regressed to maximize the log-likelihood of the marginal distribution over the sampled data [13]. For a query point a∗ the posterior prediction conditioned on the data D = {(ai , Ψ̂(ai ))}N i=1 is a normal distribution Ψ(a∗ )|D = N (µ(a∗ ), σ 2 (a∗ )), where the posterior mean µ and covariance σ 2 are computed as µ(a∗ ) = k(a∗ )K −1 Ψ̂ and σ 2 (a∗ ) = k̄(a∗ , a∗ ) − k(a∗ )K −1 k(a∗ )T , where k̄ is the kernel of the GP, k(a∗ ) = [k̄(a∗ , a1 ), ..., k̄(a∗ , aN )], Ki,j = k̄(ai , aj )+σn2 δi,j , where δi,j is the Kronecker delta, and Ψ̂ = [Ψ̂(a1 ), ..., Ψ̂(aN )]. C. Discrete-Time Lyapunov Function Consider the discrete-time dynamics (4) and a smooth reference zkref , vkref . The tracking error with respect to this reference is defined as ek = zk −zkref . We have the following standard stability result [14, Thm. 5.5]. Lemma III.2 (Discrete-Time Lyapunov Decrease). Consider the feedback vk = −Kek + vkref and the resultant feedback error dynamics ek+1 = (Ad − Bd K)ek . If there exists a function V : Rnz → R≥0 that satisfies V (0) = 0, V (ek ) > 0, V (ek+1 ) < V (ek ), (5) ∀ek ∈ Rnz \ {0}, then the closed-loop error dynamics are asymptotically stable with respect to the origin. IV. M ETHODOLOGY Our proposed architecture (Fig. 1) builds upon [11]. The control approach consists of an FMPC, a safety filter, and a system extension. The FMPC (Section IV-A) computes the optimal flat state and input to best track z ref while respecting flat-state constraints. As the nonlinear input transformation is unknown, we learn it offline from data using an individual GP for each component of v (Section IV-C). The safety filter (Section IV-D) then uses the GPs to obtain the extended system input ūk while satisfying probabilistic Lyapunov decrease, probabilistic half-space flat state constraints, and constraints on the extended input ūk and system input uk . Finally, we obtain uk from ūk in the system extension (Section IV-B). For control-affine systems, the safety filter can be written as an SOCP. As the FMPC is a quadratic program, both optimization problems are convex and can be solved efficiently. While recursive feasibility is not theoretically guaranteed, we show in our experiments that, pragmatically, feasibility is maintained. A. Flat Model Predictive Control FMPC solves a discrete-time, finite-horizon convex optimal control problem (OCP). The formulation we introduce here builds on [4] and is equivalent to the one presented in [11]. We consider the cost (zk −zkref )T Q(zk −zkref )+vkT Rvk , where Q ∈ Rnz ×nz , Q ≻ 0, R ∈ Rm×m , R ≻ 0, Q is block diagonal, and R is diagonal. The FMPC at timestep k receives the estimated flat state ẑk , has a horizon of T ∈ N and uses the first flat state zk∗ and input vk∗ of its horizon as inputs to the safety filter (Section IV-D). The OCP ensures zk∗ ∈ Z along the prediction horizon. We note that, under
the same assumptions as in [11, Sec. IV.A], Theorem 6.20 from [15] implies that for a sufficiently large T , the FMPC renders (4) asymptotically stable with respect to the reference z ref . Furthermore, from [16, Sec. 8.3], there is a closed form solution with equivalent gain matrix K ∈ Rm×n vk∗ = −K(zk − zkref ) + vkref , (6) in the absence of constraints Z. Here, K can be computed from the finite horizon discrete-time Ricatti equation (DTRE). Remark IV.1. Note that due to the block diagonal structure we impose on Q and R, we have that P from the solution of the DTRE will also be block diagonal. B. System Extension To handle flat input transformations requiring input derivatives—extending [11]—we consider a discretized dynamic extension [5], ηk = Ad,ext ηk−1 + Bd,ext ūk , uk = Cext ηk . (7) Here, ηk includes the components of the input uk and their derivatives up to the required order, ūk is the extended input, and Ad,ext and Bd,ext fulfill the exact-discretization chain integrator dynamics. Note that ūk may include both extended and non-extended input components. With this system extension (3) simplifies to vk = Ψ(zk , ūk ) = α(zk ) + β(zk )ūk . (8) for control-affine systems such as (1). This transformation Ψ is unknown and must be learned from data. Assumption IV.1. We assume that the extended input ū, the corresponding system extension (7), and the dimension of the flat state z for system (1) are known. We note that verifying the structural assumptions of our approach (flatness, relative degree, control-affine structure) requires only qualitative knowledge of the system class (e.g., from first principles or literature) and does not require detailed system identification. Remark IV.2. If (3) doesn’t require input derivatives (i.e., the system has full relative degree) then ū = u and (8) holds without system extension. C. Gaussian Process Learning We propose to learn each component of the nonlinear transformation (8) with an independent GP, as commonly done in GPMPC [2]. Thus, each of the m GPs learns: T (9) i v = i Ψ(z, ū) = αi (z) + βi (z) ū, where the index i on the bottom left denotes component i of a vector and αi and βi are the components of α and β in (8) corresponding to component i Ψ. We encode the affine structure of (9) directly in the kernel of the GP using an affine kernel, as introduced in [17], k̄(ai , aj ) = kα (zi , zj ) (10) +ūT i diag(kβ,1 (zi , zj ), ..., kβ,m (zi , zj )) ūj where ai = [zi , ūi ] is the input to the GPs. Assumption IV.2. kα and kβ,1 , ..., kβ,m are positive definite and bounded kernels. Given Assumption IV.2, the kernel (10) is also positive definite and bounded as shown in [17, Lemma 3]. Using the
affine kernel (10) the mean and covariance prediction at a query point a = [z, ū] are computed with i
µ = i γ1 (z) + i γ2 (z)T ū
i 2
i
i
T
σ = γ3 (z) + γ4 (z) ū + ū
(11) Ti
γ5 (z)ū.
(12)
i
The derivation of each γj is similar to [11], but here we consider a multi-dimensional transformation. Assumption IV.3. For each component i of the nonlinear transformation (8), the true function i Ψ(a) belongs to a Reproducing Kernel Hilbert Space (RKHS) induced by the employed kernel (10) of the corresponding GP. Furthermore, its RKHS norm is bounded as ∥ i Ψ(a) ∥k̄ ≤ Bi < ∞ for some Bi > 0. The observation noise of each GP is independent and uniformly bounded by σn . Under Assumption IV.3, [18, Thm. 6] guarantees that for some δi ∈ (0, 1), mean predictions are bounded as 1/2
Pr{∀a ∈ A, | iµ(a) − i Ψ(a)| ≤ βi i σ(a)} ≥ 1 − δi , (13) where A is compact, and β is computed based on the properties of the GP and the training data.
D. Safety Filter We combine the learned flat-input transformation with the FMPC results in a safety filter consisting of four components: probabilistic feedback linearization, a probabilistic Lyapunov decrease constraint, a probabilistic state constraint and input constraints. We then show that this safety filter is an SOCP that finds the optimal ū. To limit the search space of the optimization and to enable the stability constraint bounds, we impose ūmin ≤ ū ≤ ūmax . 1) Probabilistic Feedback Linearization: From FMPC we get the desired optimal flat states and inputs zk∗ and vk∗ . To find ūk , we minimize the distance between the expected GP prediction and vk∗ min E[∥ Ψ(zk∗ , ūk ) − vk∗ ∥2 ], ūk
where E is the expected value and Ψ = [1 Ψ, ..., m Ψ]T is the GP posterior distribution of each dimension summarized in one vector. Similar to [11], we use the definitions of the GP mean and covariance from Eqs. (11) and (12) and leave out all terms independent of ūk to get ! m X i ∗ ∗ i ∗T i ∗T 2( γ1 − i v k ) γ2 + γ4 min ūk ūk
i=1
+ūT k
m X
i
T γ2∗ i γ2∗ + i γ5∗
(14)
! ūk ,
i=1
using i γj∗ := i γj (zk∗ ). The solution to this optimization finds the ūk that, when mapped through the learned nonlinear input transformation (8), most closely matches the optimal flat input from the FMPC vk∗ . 2) Probabilistic Stability Constraint: Using the Lyapunov decrease condition (5) we formulate a probabilistic stability constraint to ensure that the input uk guarantees Lyapunov decrease of the closed loop system with learned dynamics. This will ensure that, as long as the optimization remains feasible, the tracking error will decrease. We consider a nominal controller vknom = −Kek + vkref with gain matrix K from (6) and a Lyapunov function V (ek ) = eT k P ek with P from Remark IV.1. With the nominal controller the error dynamics become
ek+1 = (Ad − Bd K) ek + Bd (Ψ(zk∗ , ūk ) − vknom ). Inserting the error dynamics into the Lyapunov decrease condition in (5) results in w1T (Ψ(zk∗ , ūk ) − vknom ) + T (Ψ(zk∗ , ūk ) − vknom ) W2 (Ψ(zk∗ , ūk ) − vknom ) ≤ w3 which can be written as T (Ψ(zk∗ , ūk ) − vknom + w4 ) W2 (Ψ(zk∗ , ūk ) − vknom + w4 ) 1 ≤ w3 + w1T W2−1 w1 , (15) 4 T T where w1T = 2eT k (Ad − Bd K) P Bd , W2 = B d P Bd , T T w3 = ek P − (Ad − Bd K) P (Ad − Bd K) ek − ε, w4 = 12 W2−1 w1 , and ε > 0 is a small constant to allow the inequality to be non-strict. Theorem IV.1. Given (13) and Remark IV.1 the Lyapunov decrease condition (15) holds when m X 2 1/2 i ∗ nom + i w4 + Li βi i σ(zk∗ , ūk ) i W2 µ(zk , ūk ) − i v k i=1
1 (16) ≤ w3 + w1T W2−1 w1 , 4 Pm i with probability δs = 1 − i=1 δi . Here iµ and σ 2 are the GP mean and covariance prediction of i Ψ, i W2 is the diagonal component i of W2 and Li is given by Li = min max max 2 i W2 (| iµ(zk∗ , s) − i v nom + i w4 | + ρ̄ i σ(zk∗ , s)), k s∈[ūk
,ūk
]
where ρ̄ is the Gaussian quantile function at probability δ̄. Proof: From Remark IV.1 P is a positive definite blockdiagonal matrix. Combined with the block structure of Bd in (4) W2 is a diagonal matrix W2 = diag(1 W2 , ..., m W2 ). Thus we can write the left hand side of (15) component wise m X 1 2 nom + i w4 ) ≤ w3 + w1T W2−1 w1 , i W2 (i Ψ − i v k {z } | 4 i=1 i g(i Ψ)
where i g(i Ψ) is quadratic in i Ψ (arguments omitted for readability). Since ūk is bounded, i g(i Ψ) is also bounded and i g is Lipschitz continuous with constant L̃i = + i w4 ) |, which we bound by Li maxi Ψ |2 i W2 (i Ψ − i v nom k with high probability δ̄. Applying (13) gives | i Ψ − iµ| ≤ 1/2 βi i σ with probability 1P − δi and Boole’s inequality yields m the joint probability 1 − i=1 δi . We use the definitions of the GP mean and covariance Eqs. (11) and (12) to write (16) solely as a function of ūk q m X 1/2 i ∗ T i ∗ γ3 + i γ4∗ ūk + ūT Li βi k γ5 ūk ! i=1 m X i ∗i ∗ T + ūT ūk (17) k i W 2 γ 2 γ2 ! i=1 m X T − i W2 2i γ1∗ + 2(i w4 − i v nom ) i γ2∗ ūk ≤ k i=1 m X 2 1 i ∗ nom + w3 + w1T W2−1 w1 − + i w4 i W2 γ1 − i v k 4 i=1 3) Probabilistic State Constraint: Because the safety filter can modify the flat input, and the learned mapping Ψ(zk∗ , ūk ) is uncertain, we need a constraint in the safety filter to ensure zk+1 ∈ Z with high probability. With the GP predictions and the discrete-time linearized dynamics (4) the mean of the next state becomes
µz,k+1 = Ad zk∗ + Bd [1µ(zk∗ , ūk ), ..., mµ(zk∗ , ūk )]T . For brevity, we only consider uncertainty in µz,k+1 due to i 2 σ , i ∈ {1, ..., m}. Therefore the uncertainty of the next state 1 2 ∗ m is Σz,k+1 = Bd diag σ (zk , ūk ), ..., σ 2 (zk∗ , ūk ) BdT . The off-diagonal terms in the covariance matrix are zero because the components of the flat input vector are assumed to be independent of each other. Using Σz,k+1 we tighten Z using probabilistic reachable sets (PRS). Following the one-step PRS tightening from [11] and [2], the tightened constraint set Z R (Σz,k+1 ) := {zk+1 |hT j zk+1 ≤ bj − ρ(1 − q
δj ) hT j Σz,k+1 hj }, for j = 1, . . . , nc , ensures that under R the uncertain Pncdynamics µz,k+1 ∈ Z with probability level δc = 1 − j=1 δj . Note ρ(1 − δj ) is the quantile function of a standard Gaussian random variable at the probability 1−δj for some small δj . 4) Input Constraints: In addition to the constraints on the extended input ūmin ≤ ū ≤ ūmax we consider constraints on the system input u. This is necessary to account for actuator limits in practice. To constrain u we use the system extension (7) to relate ū to the bounds on u as umin ≤ Cext (Ad,ext ηk−1 + Bd,ext ūk ) ≤ umax (18) 5) Safety Filter as a SOCP: Here we show that the safety filter takes the form of a SOCP. Theorem IV.2. The safety filter consisting of probabilistic feedback linearization (14), a probabilistic Lyapunov decrease constraint (17), probabilistic state constraint (Section IV-D.3) and input constraints (18) is written as a SOCP with " ! # m X i ∗ ∗ i ∗T i ∗T min 2( γ1 − i v k ) γ2 + γ4 , 1, 0 ũ ũ
i=1
s.t. ∥ Āi ũ − b̄i ∥≤ c̄i ũ + d¯i i ∈ {1, 2, 3, 4, . . . , nc }, umin ≤ Cext (Ad,ext ηk−1 + Bd,ext ūk ) ≤ umax , ūmin ≤ ū ≤ ūmax ,
T where ũ = [ūT k , q1 , q2 ] , and q1 and q2 are auxiliary variables. The optimal input ūk is within the input bounds and guarantees Lyapunov function decrease and zk ∈ Z with probability of at least δs + δc − 1. Proof: First, we introduce q1 to be able to write the cost function as a linear function in ũ. Take q1 ≥ Pm i ∗ i ∗ T i ∗ ūT γ γ + γ ū 2 2 5 k . Using the same steps as in k i=1 our previous work [11] we canreformulate this as (1+q1 )2 ≥ P m i ∗ i ∗T 4ūT + i γ5∗ ūk + (1 − q1 )2 which can be k i=1 γ2 γ2 written as a standard SOC constraint using the Cholesky decomposition of i γ5∗ = i L∗ i L∗T to give expressions for Ā1 , b̄1 , c̄1 and d¯1 . Next, we convert the stability constraint (17) as two SOC constraints using q2 . We write Pm i ∗i ∗ T ūT W γ γ ū ≤ q2 as a SOC constraint 2 k k i=1 i 2 2 T with b̄2 = [0, ..., 0, 1] , c̄2 = [0T , 0, 1], d¯2 = 1, and Ā2 is of zeros the first m columns as p with p a matrix 1 ∗T m ∗T T col W γ , . . . , W γ , 0 , and whose last 2 2 1 2 m 2 entry is −1. Using the subadditivity of the square root, the SOC constraint for (17) with Ā3 as a zero matrix with the first few columns as
T 1/2 1/2 T col L1 β1 1 L∗ , 0 , . . . , Lm βm [m L∗ , 0] , 0 and the last entry as 1; with b̄3 = col(b̄3,1 , . . . , b̄3,m , 0) where 1 i ∗ −1 i ∗ γ4 q −2 L i γ3∗ 1 i ∗ −1 i ∗ 2 γ4 ) 1/2 m − 1( 2 L = Li β 1/2 b˜i b̄3,i = Li βi i .. . q iγ ∗ 3 1 i ∗ −1 i ∗ 2 γ4 ) m − m( 2 L hP i i ∗T m i ∗ nom c̄3 = − W 0 0 2 γ + 2( w − v ) γ 2 1 4 2 i i i k i=1 m X 2 1 T −1 i ∗ nom ¯ d3 = w3 + w1 W2 w1 − + i w4 i W2 γ1 − i v k 4 i=1
where i (·)2 is the square of component i of (·). We bound the terms in b̄3,i such that they remain real. For the probabilistic state constraint in Section IV-D.3 we split up hj into components of the same size as the block vectors of B qdPsuch that the constraint becomes m 2 i 2 ∗ hT µ ≤ b − j j z,k+1 i=1 ws,i σ (zk , ūk ) where ws,i = q ρ(δ) i hj T bd,i bT d,i i hj . This has the same form as the probabilistic stability constraint without auxiliary variables and can therefore be written as a SOC constraint using the same procedure as in [11], to give expressions for the SOC constraints i = 4, . . . , nc . Applying Boole’s inequality yields a joint satisfaction probability of at least δs + δc − 1 for both the Lyapunov decrease and flat-state constraints when applying ūk to (1). Remark IV.3. Because the Lyapunov function decreases with high probability between successive sampling instants, if the problem remains feasible then ek → 0 as k → ∞. Although our formulation doesn’t theoretically guarantee recursive feasibility, we demonstrate in Section V that, practically, feasibility is often maintained. V. S IMULATION AND H ARDWARE E XPERIMENTS We evaluate our control approach in both simulation and hardware on an unconstrained and state/input-constrained lemniscate trajectory tracking task. To highlight our performance relative to idealized dynamics, we compare our approach with perfect-knowledge FMPC and NMPC; we also compare to GPMPC [2], [19] as a marker of stateof-the-art learning-based control. All code is available at github.com/utiasDSL/mimo fmpc socp. The 2D quadrotor model has state x = [x, ẋ, z, ż, θ, θ̇]T and input u = [Tc , θc ]T . Here (x, y) is the position, θ is the roll angle, Tc is the thrust, and θc is the commanded roll angle. The system dynamics are ẍ = sin θ(β2 + β1 Tc ), z̈ = cos θ(β2 + β1 Tc ) − g and θ̈ = α1 θ + α2 θ̇ + α3 θc . We used a Bitcraze Crazyflie 2.1 as our hardware platform; its identified model parameters can be found in our code. This system is flat with flat output y = (x, z)T , and the flat state vector is z = [x, ẋ, ẍ, x(3) , z, ż, z̈, z (3) ]T . As the system has non-full relative degree, we extend the system with a double integrator on Tc and the extended input becomes ū = [T¨c , θc ]T . A. Simulation Experiment and Results All controllers were tuned to maximize performance while maintaining reasonable computational efficiency. All controllers are run at 100 Hz with a horizon of 50 using ACADOS [20], a state-of-the-art solver for embedded optimal control,
ensuring a fair comparison across formulations. The gains Q and R for each method were kept as similar as possible while optimizing individual performance. All GPs used squared exponential kernels and were trained offline on 600 data points. For GPMPC, two GPs were used to model error residuals to account for incorrect prior parameters. The first GP corrects β2 + β1 Tc and the second corrects θ̈. To keep GPMPC tractable, sparse GP inference with 75 inducing points was used. For our FMPC+SOCP method, the points {zk , ūk , vk } were sampled around the reference. In both the constrained and unconstrained tasks, the simulation was run 30 times starting at a randomized state within a small ball around the reference. On held-out data, all GPs achieved relative RMSE below 2% of the output range, with over 96% of test points falling within the 2-σ confidence bounds, except for the first flat input component (86%), which is stiff and more difficult to calibrate. In Fig. 2a, we show a comparison of the tracking error without constraints. After a transient phase, our approach with purely learned dynamics achieves tracking performance similar to FMPC and slightly better than NMPC and GPMPC. The observed transient error arises from the limits on ū, which restricts rapid thrust changes. These are even more significant when starting from hover, as seen in the hardware comparison in Fig. 2b. In Fig. 2c, we demonstrate the behavior under constraints on the flat state x and on the input Tc . Our approach anticipates the state constraint on the FMPC horizon and thus satisfies it with minimal tracking error. In contrast, the input constraint on Tc is enforced in the single-step safety filter (Fig. 3), so it cannot be anticipated. This can cause a slight 2 overshoot bounded by ∆Tc = (Ṫc,0 )/(2T̈c,max ), where Ṫc,0 is the thrust derivative when the constraint becomes active and T̈c,max is the bound on the corresponding component of ūk ; we tighten the constraint by this amount. Minor residual violations can also arise from ACADOS solver instabilities, as observed for both NMPC and our method. The distribution of inference times on the constrained task of each approach are shown in Fig. 4. We see that our method averages 3.1 ms, which is less than 3.4 ms for NMPC and more than three and a half times faster than 11.2 ms for GPMPC. The larger variance in NMPC and GPMPC comes from the nonlinear optimizations that sometimes— particularly near constraint boundaries, far away from warmstart solutions, or in uncertain regions—require significant compute. These results confirm that our learning-based flat control delivers significant computational benefits without sacrificing performance, even when accounting for dynamics uncertainty and constraints. B. Hardware Experiments and Results Hardware experiments were performed using a Bitcraze Crazyflie 2.1 quadrotor. A Vicon motion capture system with an extended Kalman filter was used to estimate the system state x and flat state z assuming perfect system knowledge. The GPs trained in simulation were used for the hardware experiments—training the GPs on real-world data is left for future work. The controllers were run at 50 Hz on a standard laptop from a base station and telemetry was wirelessly sent to the drone. Instead of starting near the reference as in the simulation tasks, the drone started from hover.
1.6
FMPC+SOCP (0.007 m) GPMPC (0.018 m)
0.03 0.02 0.01 0.00
FMPC+SOCP (0.117 m) 0.3
0.2
0.1
5
10
15
20
1.4
1.2 1.0 0.8 0.6
1.0 0.8 0.6
0.4 0
25
1.2
FMPC+SOCP
0.0 0
1.6
1.4
Position z (m)
0.04
FMPC (0.095 m) Position z (m)
FMPC (0.007 m)
NMPC FMPC+SOCP GPMPC
NMPC (0.076 m)
0.4
NMPC (0.006 m)
Tracking error (m)
Tracking error (m)
0.05
2
4
6
8
10
12
14
Time (s)
Time (s)
(a) Simulation
(b) Hardware
0.4 −1.0
−0.5
0.0
0.5
1.0
−1.0
−0.5
Position x (m)
(c) Simulation
0.0
0.5
1.0
Position x (m)
(d) Hardware
Angle θc (rad)
Thrust Tc (N)
Fig. 2. Comparisons of the position tracking errors on the unconstrained flat state tracking task for the (a) simulated and (b) hardware show that our formulation is very competitive with the controllers using the true dynamics. When including constraints, we can see the performance is still competitive in (c) simulation, and that on (d) hardware our controller still successfully respects the constraint. All simulation plots show the mean (solid) and 95% confidence interval (shaded) across all the initial conditions. The red shaded region represents the flat state constraint. The root-mean-square error of each controller in the unconstrained case is given in the legend inside the parenthesis. 0.4
on learning-based FMPC with theoretical recursive feasibility and asymptotic stability under constraints.
0.3 0.2
R EFERENCES
0.1 2
0
5
10
0
5
10
0
15 NMPC FMPC+SOCP GPMPC
20
25
20
25
−2 15
Time (s)
Inference Time (ms)
Fig. 3. The mean input (solid) and 95% confidence interval (shaded) across all the initial conditions for the constrained simulation task. The red shaded region represents the constraint. 12 10
Mean Median
8 6 4 2 NMPC
GPMPC
FMPC+SOCP
Fig. 4. Violin plot of the average inference times on the constrained task. Outliers were filtered using 1.5 times the interquartile range.
For the unconstrained experiment, the tracking error in Fig. 2b indicates that our approach is competitive with the NMPC and FMPC. On the constrained tracking task, our proposed approach respects all constraint boundaries and remains feasible for all time (Fig. 2d). Due to model mismatch, state-estimation noise, and the use of GPs trained solely in simulation, the tracking errors are larger than in simulation for all approaches. Offset-free MPC techniques [21] could further reduce steady-state bias in hardware and are left for future work. VI. C ONCLUSION We introduced a computationally efficient learning-based control method for nonlinear multi-input, control-affine differentially flat systems that guarantees probabilistic Lyapunov decrease and probabilistic half-space flat state constraint satisfaction and respects input constraints. While recursive feasibility is not yet theoretically guaranteed, we show through simulation and hardware experimentation that our controller remains practically feasible, stable, and respects constraints. Scalability is aided by the convexity and block-diagonal decoupling of both subproblems, and perdimension GP inference. Overall, our experimental results show that our flatness-exploiting learning-based controller is significantly more efficient than GPMPC, the gold-standard learning-based controller, and doesn’t compromise on performance. This work lays a strong foundation for future work
[1] L. Brunke, et al., “Safe Learning in Robotics: From Learning-Based Control to Safe Reinforcement Learning,” Annu. Rev. Control Robot. Auton. Syst., vol. 5, no. 1, pp. 411–444, 2022. [2] L. Hewing, J. Kabzan, and M. N. Zeilinger, “Cautious Model Predictive Control Using Gaussian Process Regression,” IEEE Trans. Contr. Syst. Technol., vol. 28, no. 6, pp. 2736–2743, 2020. [3] M. Fliess, et al., “Flatness and Defect of Non-Linear Systems: Introductory Theory and Examples,” Int. J. Control, vol. 61, no. 6, pp. 1327–1361, 1995. [4] M. Greeff and A. P. Schoellig, “Flatness-Based Model Predictive Control for Quadrotor Trajectory Tracking,” in IEEE Int. Conf. Intell. Robots Syst. (IROS), 2018, pp. 6740–6745. [5] A. Isidori, Nonlinear Control Systems. London: Springer, 1995. [6] S. Sun, et al., “A Comparative Study of Nonlinear MPC and Differential-Flatness-Based Control for Quadrotor Agile Flight,” IEEE Trans. Robot., vol. 38, no. 6, pp. 3357–3373, 2022. [7] M. Greeff and A. P. Schoellig, “Exploiting Differential Flatness for Robust Learning-Based Tracking Control Using Gaussian Processes,” IEEE Control Syst. Lett., vol. 5, no. 4, pp. 1121–1126, 2021. [8] B. Akbari and M. Greeff, “A Computationally Efficient LearningBased Model Predictive Control for Multirotors under Aerodynamic Disturbances,” in Int. Conf. Unmanned Aircr. Syst. ICUAS, 2024. [9] B. Akbari, J. Frank, and M. Greeff, “Tiny learning-based mpc for multirotors: Solver-aware learning for efficient embedded predictive control,” Mechatronics, vol. 115, p. 103452, 2026. [10] M. Greeff, A. W. Hall, and A. P. Schoellig, “Learning a Stability Filter for Uncertain Differentially Flat Systems using Gaussian Processes,” in Proc. IEEE Conf. Decis. Control (CDC), 2021, pp. 789–794. [11] A. W. Hall, M. Greeff, and A. P. Schoellig, “Differentially Flat Learning-Based Model Predictive Control Using a Stability, State, and Input Constraining Safety Filter,” IEEE Control Syst. Lett., vol. 7, pp. 2191–2196, 2023. [12] J. Levine, Analysis and Control of Nonlinear Systems: A Flatnessbased Approach. Berlin, Heidelberg: Springer, 2009. [13] C. E. Rasmussen and C. K. I. Williams, Gaussian Processes for Machine Learning. Cambridge, MA, USA: MIT Press, 2006. [14] K. Ogata, Discrete-Time Control Systems, 2nd ed. Upper Saddle River, NJ, USA: Prentice-Hall, 1998. [15] L. Grüne and J. Pannek, Nonlinear Model Predictive Control. Cham: Springer International Publishing, 2017. [16] F. Borrelli, A. Bemporad, and M. Morari, Predictive Control for Linear and Hybrid Systems. Cambridge: Cambridge University Press, 2017. [17] F. Castañeda, et al., “Gaussian Process-based Min-norm Stabilizing Controller for Control-Affine Systems with Uncertain Input Effects and Dynamics,” in Proc. Am. Control Conf., 2021, pp. 3683–3690. [18] N. Srinivas, et al., “Information-Theoretic Regret Bounds for Gaussian Process Optimization in the Bandit Setting,” IEEE Trans. Inform. Theory, vol. 58, no. 5, pp. 3250–3265, 2012. [19] Z. Yuan, et al., “Safe-Control-Gym: A Unified Benchmark Suite for Safe Learning-Based Control and Reinforcement Learning in Robotics,” IEEE Robot. Autom. Lett, vol. 7, no. 4, pp. 11 142–11 149, 2022. [20] R. Verschueren, et al., “acados – a modular open-source framework for fast embedded optimal control,” Math. Program. Comput., 2021. [21] M. Morari and U. Maeder, “Nonlinear offset-free model predictive control,” Automatica, vol. 48, no. 9, pp. 2059–2067, 2012.