The MRS UAV System: Pushing the Frontiers of Reproducible Research, Real-world Deployment, and Education with Autonomous Unmanned Aerial Vehicles

In plain words

A small flying robot with four or more propellers is unstable by nature, so it needs a control loop that corrects itself many times per second. This paper describes a complete software stack for such drones, released as open source, and built by a robotics group that used it in real competitions. The authors wanted a system that can be copied, rerun and trusted outside a lab.

The stack has three main parts. First, a set of simple position filters runs at the same time, each fed by a different mix of sensors such as GPS, a laser or camera odometry, and a heading estimate. The drone can switch between them in flight, so it can move from GPS to an indoor source without jumping. Second, the controller turns the desired path into a force to produce, and then into motor commands. It also learns the constant pushes it feels, such as wind or a battery whose mass changed, and cancels them. Third, a second controller uses a short-horizon optimisation to stay smooth when the position estimate is noisy. The authors describe the heading of the drone as a direction vector rather than with yaw angles, which they found easier to reason about.

TL;DR

Open-source, ROS-based full-stack multirotor UAV system from CTU Prague. It has a bank of linear Kalman filters for several frames of reference, a heading-based SO(3) force tracker, and a linear MPC force feedback with disturbance integrators. The results are qualitative (simulation vs. real flight on a Tarot 650) plus circle-tracking errors and a jump-handling test. There is no benchmark table.

Problem

Research on multirotor control is often tested only in simulation, or on a platform whose software cannot be transferred to other hardware §1.1. Many open-source UAV stacks depend on one flight controller's built-in position control, or on a GPS-only setup, so they do not work outside a lab §1.1. The authors also argue that controllers which assume a clean state estimate can become unstable when the estimate is noisy or jumps, and that such jumps can push the drone into unplanned motion §3.5. The paper treats these as the practical problems for a deployed system §1.

Method

The system has four parts §1.3, fig. 3: a mission layer that supplies a position and heading goal; a reference tracker that turns it into a smooth, feasible, full-state reference at 100 Hz; a feedback controller that outputs an attitude rate and a collective thrust; and a state estimator that runs a bank of filters. The embedded Pixhawk flight controller (PX4) closes the attitude-rate loop §1.3, §3. The drone model is a rigid body with translation governed by thrust and gravity and rotation given by a rotation matrix, with the heading defined from the first body axis §2.

Heading convention. Heading is the azimuth of the first body axis b1 in the world plane, η = atan2(b1·e2, b1·e1). The authors advise against Tait-Bryan or Euler yaw because the yaw depends on the other two rotations §2. The translational estimator, reference tracker and controllers all carry this heading η and its rate.

State estimation (§3). Each filter is a discrete, decoupled linear model with 11 states: position, velocity and acceleration for x, y and z, plus heading and heading rate §3.1, state vector. The model has three first-order transfer parameters p1 (horizontal), p2 (vertical) and p3 (heading), with the two horizontal axes assumed identical §3.1. The filters are standard linear Kalman filters, and the authors choose them over EKF/UKF for cost and tuning reasons §3.2. The heading rate is taken from the total derivative of atan2 applied to the rotation matrix, not from the yaw rate alone; the authors state this is only valid when the drone is not nearly horizontal §3.4. After estimation the orientation is re-based by rotating the Pixhawk rotation matrix about e3 by the difference between the estimated and original headings §3.3.

Bank of filters (§3.5). Several filters run in parallel, each corrected by a subset of the sensors. One filter is the primary hypothesis used for feedback. An arbiter switches the primary when requested or when the current filter is unreliable, and otherwise chooses the filter with the smallest trace of its covariance §3.5. Each switch changes the coordinates of the drone without moving it. The controllers and trackers are then transformed to the new frame §3.5, §7.2.

Reference tracking (§4). The MPC tracker is a linear MPC on a virtual point-mass model. It outputs a full-state reference at 100 Hz: position and its derivatives up to jerk, heading and heading rate §4.1. Take-off uses an admittance rule that limits how far the reference can move from the current state, and landing sets the reference altitude below the estimated altitude §4.2. A speed tracker passes through a desired velocity or acceleration for swarming §4.3.

Feedback control (§5). The base is the force-tracking part of the geometric controller by Lee et al. §5.1, §1.1. The desired force is built from position, velocity and feedforward terms, with the desired heading vector fixing the orientation §5.1, §5.6. Two ways of building the desired rotation matrix are given: the original Lee construction R_d^o, and a heading-compliant construction R_d based on an oblique projection of the desired heading vector onto the plane orthogonal to the desired thrust direction §5.1.1, §5.1.2. Rotation error uses the standard antisymmetric part of R_d^T R − R^T R_d §5.1.2. The desired angular rate is −k_R ∘ e_R + ω_j − ω_c, where ω_j is a feedforward term from the desired jerk and ω_c compensates the heading rate caused by roll and pitch rates §5.1.2, §5.2, §5.3. The yaw-rate feedforward is recovered from the desired heading rate by a projection §5.2. Collective thrust follows a quadratic model, T_d = a_t·sqrt(f_d) + b_t, fitted by least squares on measured thrust and mass pairs §5.4.

Disturbance estimation (§5.5). Two integrators estimate force disturbances without using the drone model or the state estimator. One integrates the position error in the world frame. The other integrates the error in the body frame, using the heading at each sample, and rotates the sum by the current heading §5.5. The paper splits the disturbance equally between the two frames until the heading changes §5.5. The estimated mass is the nominal mass plus the vertical component of the total disturbance §5.5. The output is used for both feedback controllers and for an "unbiased" desired acceleration §5.8.

SE(3) controller (§5.6). The agile controller adds position and velocity feedback, reference acceleration, gravity and disturbance compensation to the desired force §5.6. The horizontal disturbance terms are applied, and the vertical one is absorbed in the estimated mass §5.6.

MPC force feedback (§5.7). A linear MPC runs on a nine-state model per axis, with acceleration as the control input and no input penalty [§5.7.1, §5.7.2]. The model step is 0.05 s §5.7.1. The cost sums squared tracking errors over a horizon of length n with weights Q and S, and constrains state bounds and input slew rate. The slew rate constrains the jerk §5.7.2. The first input is the desired acceleration c_d §5.7.2. The desired force adds the reference acceleration, c_d, gravity and disturbance compensation §5.7. The controller follows the same quadratic-programming formulation at 100 Hz §5.7.2.

# MPC force feedback, one control step  [§5.7] [§5.7.2]
# Input:  estimate state, reference chi, disturbances, mass estimate
# Output: desired angular rate (3,), desired thrust force (scalar)

def mpc_feedback(state, reference, d_w, d_b, mass_e, R, eta_d):
    # 1. Initial condition: estimate, or reference if estimate violates limits  [§5.7.2]
    x_m0 = initial_state(state, reference)
    # 2. Solve the QP over the horizon; keep the first acceleration  [§5.7.2]
    c_d = solve_mpc(x_m0, reference)
    # 3. Desired force: reference + MPC acceleration + gravity - disturbances  [§5.7]
    f_d = (mass_e * (reference.acc + c_d) + mass_e * g * e3
           - (d_w + d_b) * mask_xy)                    # mask_xy = [1, 1, 0]
    # 4. Desired rotation that keeps the heading eta_d  [§5.1.2]
    R_d = heading_compliant_rotation(f_d, eta_d)
    # 5. Rotation error  [§5.1.2]
    e_R = rotation_error(R_d, R)
    # 6. Feedforward rate from desired jerk, and parasitic heading-rate term  [§5.1.2] [§5.3]
    omega_j = feedforward_rate(f_d, mass_e, R_d, reference.jerk)
    omega_c = parasitic_rate(R, omega_d_prev)
    omega_d = -k_R * e_R + omega_j - omega_c   # gains: not stated
    # 7. Scalar thrust along the body axis  [§5.1.2]
    f_scalar = dot(f_d, R[:, 2])
    return omega_d, f_scalar

# Disturbance integrators, one update per control step  [§5.5]
# Input:  position error e_p (3,), heading eta, time step dt
# Output: world and body disturbances d_w (3,), d_b (3,)

def update_disturbances(e_p, eta, dt, acc):
    # 1. World-frame integral of the position error  [§5.5]
    acc.world = acc.world + k_iw * e_p * dt          # elementwise
    # 2. Body-frame integral: error rotated into the body frame at this sample  [§5.5]
    H_n = heading_rotation(eta)                      # yaw-only rotation
    acc.body = acc.body + k_ib * (H_n.T @ e_p) * dt  # elementwise
    # 3. Rotate the accumulated body sum by the current heading  [§5.5]
    d_w = acc.world
    d_b = heading_rotation(eta) @ acc.body
    return d_w, d_b

Unbiased acceleration (§5.8). The controller reports the acceleration it produces after removing gravity, the body and world disturbances, and the mass difference. This keeps the state estimator's model free of biases [§3.1.1, §5.8].

Failsafe (§5.10). If position feedback is lost, a feedforward controller keeps the drone level and reduces thrust from the hover-thrust value at a fixed rate, using the last estimated mass §5.10.

Trajectory generation (§6). Waypoint paths are turned into polynomial trajectories with a nonlinear time-optimisation method that extends a prior approach; only segments that violate limits are stretched §6.

Hyperparameters

Results

Author-reported. Setup: a Tarot 650 multirotor, flown in the real world and in the authors' Gazebo simulator; the simulator is claimed to match the real flights closely §8. No benchmark table is given. Quantitative values are in the text and captions:

Other sections (competitions, swarming, mutual UAV detection) are narrative examples and carry no comparable numbers §9.

Paper vs code

Checked: the official wrapper repository ctu-mrs/mrs_uav_system (linked in the paper) at commit ef41189e981fcfdbd7512db93f0bf6f936e2fa7b, its default branch, which is the ROS 2 line. It contains only build files and documentation, under BSD-3-Clause LICENSE. The controller code is in ctu-mrs/mrs_uav_controllers, which the wrapper does not contain. For that repository I read the ros1 branch via WebFetch, not a pinned clone, so its commit is unknown. The files read were config/public/se3_controller.yaml, config/public/mpc_controller.yaml and src/mpc_controller.cpp. Line numbers from those files were not verified, so they are not cited. The MPC solver is distributed as a prebuilt library, lib/MpcControllerSolver/x64/, and its source was not available to me. Not checked: the bank-of-filters and Kalman code (in mrs_uav_odometry), the thrust model values, the failsafe, and trajectory generation.

• Rotation option: paper states it prefers the original Lee construction R_d^o §5.1.2. The official SE(3) config defaults rotation_matrix to 1, the "baca" oblique projection, i.e. R_d se3_controller.yaml.
• MPC speed limit: paper says the controller is limited to 2 m/s, 2 m/s² and 5 m/s³ §5.7.2. The code passes 1.5 × the constraint group's horizontal speed to the solver mpc_controller.cpp, so the effective bound depends on the active constraint group §7.5.
• Initial state: paper uses the estimate if constraints hold, and otherwise the reference's velocity and acceleration §5.7.2. The code always takes position from odometry and switches acceleration and velocity to the reference only when they exceed 1.5 × the limit mpc_controller.cpp.
• Penalties: the paper gives Q and S as values §5.7.2. The code loads Q and S from parameters under mpc/mpc_parameters/horizontal and vertical, and zeroes the position entries when position feedback is off mpc_controller.cpp. The values were not found in the public YAML I read.
• Time step: the paper's model step is 0.05 s §5.7.1 while the loop runs at 100 Hz §5.7.2. The code takes the loop dt from message timestamps and uses dt1/dt2 from parameters, whose values were not found [code: src/mpc_controller.cpp].
• Solver literal: the code passes a literal 999 as a second argument to the solver limits, meaning not documented in the code I read mpc_controller.cpp.
• Mass estimator and integral gains: the paper states no values §5.5. The ros1 public config sets kiw = kib = 0.1 with limits 10, and a mass estimator rate km = 0.5 kg/s with a limit of 20 kg mpc_controller.yaml.
• Parasitic heading compensation: the paper describes it §5.3. The SE(3) config enables parasitic_pitch_roll se3_controller.yaml.
• Version: the paper describes the 2020 ROS 1 system. The wrapper's default branch is a ROS 2 port, so the code read here is not necessarily the paper's exact version.

Pitfalls

• Two time steps. The MPC model uses Δt = 0.05 s §5.7.1 while control runs at 100 Hz §5.7.2. Do not reuse one step for both. Using the 0.05 s model step in the 100 Hz loop, or the reverse, changes the predicted states.
• Units of the disturbance terms. The paper gives disturbances in newtons and the unbiased acceleration in m/s² §5.5, §5.8. Divide by the estimated mass before adding the acceleration.
• Which rotation matrix. The paper states a preference for R_d^o §5.1.2 but its own heading-compliant R_d is the one that keeps heading during manoeuvres. The official SE(3) default is the oblique construction. Pick one and test it.
• Heading singularity. The heading and its rate are defined only when the first body axis is not near the vertical, |e3·b1| < 1 §2, §3.4. Clamp or switch before the tilt gets large.
• Yaw rate versus heading rate. Using the yaw rate as the heading rate is only an approximation for tilts up to about 10° §3.4.
• Integral gains and feedback gains are not given. k_R, k_p, k_v, k_iw and k_ib are not stated in the paper §5. The ROS 1 config values are one source, but they are not the paper's.
• Estimator frame switches. An arbiter switch moves the reported coordinates without moving the drone, so controllers must transform their internal state §3.5, §7.2. Skipping this produces a jump.
• Mass in the MPC. The MPC desired force uses the estimated mass m_e, so a stale mass estimate biases the output §5.7.
• Figure-based claims. Several results (noise, jump, sim vs. real) are only shown in figures, and the paper gives no metric or variance for them fig. 12 fig. 13.
• Revision markup. The arXiv v6 source keeps \added/\updated revision macros. This card uses the wording as it reads without them.

Lineage

• Builds on: SE(3) geometric tracking, Lee et al. 2010 (cited as [lee2010geometric] in the paper) for the force-tracking structure §5.1. The MPC tracker is from Baca et al. 2018 (cited as [baca2018model]) §4.1. The trajectory method extends Richter et al. 2016 (cited as [richter2016polynomial]) §6.
• Compared against: RotorS, OpenUAV, XTDrone, Aerostack, ReCOPTER and MAVwork in the paper's comparison tables tab. 2 and in §1.1.
• Superseded by: not known from the wiki.