Crazyflie hover control: LQR, MPC and DeePC

Quadrotor simulation for a Crazyflie-sized vehicle, linearized about hover at [0, 0, 1], with three controllers: an LQR, a model-based MPC, and a data-driven predictive controller (DeePC) that uses only recorded input/output data. The attitude is a unit quaternion, so the 13-state dynamics are reduced to a 12-state error representation through the attitude Jacobian.

Layout

src/                the modules below
data/               Hankel matrices built by collect_data.py (regenerable, gitignored)
images/             generated plots
  setpoint_0.5/       DeePC flying to the setpoint [0.5, 0.5, 0.5]
  regulate_hover/     DeePC regulating back to hover from a circle of start points
requirements.txt    dependency list

Files

file	what it is
src/quaternion_utils.py	quaternion helpers in jax.numpy: hat, L, qtoQ, G, rptoq, qtorp, attitude Jacobian E; also enables jax x64
src/quadrotor.py	Quadrotor class — 13-state dynamics, jitted RK4 step, jacobians() by JAX autodiff
src/quadrotor_lqr.py	linearizes about hover, solves the Riccati equation for K, defines the shared error_state, simulates with thrust clamping
src/quadrotor_qp.py	MPC: the linear model in a horizon-20 QP with per-motor thrust limits (CVXPY/Clarabel)
src/collect_data.py	builds and compares two Hankel matrices, saves them to data/
src/quadrotor_deepc.py	DeePC to a setpoint, loading the saved matrices
src/deepc_regulate.py	DeePC regulating to hover from several initial conditions, both datasets
src/plot.py	all figures; resolves images/ from its own file location, not the working directory

Setup

python -m venv .venv
source .venv/bin/activate
pip install -r requirements.txt

Dependencies are numpy, scipy, jax, cvxpy, clarabel and matplotlib. The CPU build of jaxlib is enough — nothing here needs a GPU.

The environment is activated per shell — run source .venv/bin/activate again in each new terminal, and deactivate to leave it. In VS Code, also run Python: Select Interpreter and choose the one under .venv; the Interactive Window uses the selected interpreter rather than whatever your terminal has activated, and mismatching the two produces confusing import errors.

Running

The model-based controllers stand alone:

python src/quadrotor_lqr.py    # writes images/quadrotor_lqr.png
python src/quadrotor_qp.py     # writes images/quadrotor_qp_mpc.png

DeePC needs the Hankel matrices built first. collect_data.py writes data/*.npz; the controllers load it rather than importing collect_data, so they do not re-run the rollouts on every invocation:

python src/collect_data.py                                    # once; writes data/ + images/hankel_comparison.png
python src/quadrotor_deepc.py                                 # setpoint task, gaussian dataset
python src/quadrotor_deepc.py --dataset lqr --stride 8 --filter-dead 0.5
python src/deepc_regulate.py                                  # regulation task, both datasets

quadrotor_deepc.py options:

flag	meaning
--dataset {rnd,lqr}	rnd = gaussian excitation at hover, lqr = LQR trajectories from many initial states
--stride N	keep every Nth Hankel column. Fewer columns means a faster QP
--filter-dead FRAC	keep only the most active FRAC of columns before striding (this is for the LQR with multiple initial conditions dataset. Filter out redundant states at hover)

deepc_regulate.py takes --radius, --n-ic and --seconds.

Any script can be run from the repository root or from inside src/.

The two datasets, and which task suits them

collect_data.py builds two Hankel matrices, both from closed-loop data with the hover LQR running:

• rnd — the vehicle held at hover with gaussian noise added to the motors: a small independent per-motor term (torque) plus a common-mode term (net thrust).
• lqr — 40 trajectories from displaced initial states converging back to hover, with a small amount of Gaussian noise added per motor.

Which one is better depends entirely on the task, and the ordering reverses:

dataset	fly to [0.5, 0.5, 0.5]	regulate to hover
rnd	0.0087 m	0.0051 m
lqr	0.0750 m	0.0020 m

Regulating to hover is the task the data actually describes — the lqr set consists of trajectories doing exactly that, so it is interpolation. The setpoint is 0.866 m from hover while the furthest recorded window ends 0.075 m away, so that task is extrapolation to a state no recorded window ever visits, and the dataset collected at hover does better there.

Notes

rank() is necessary but nowhere near sufficient. Across every sweep in this repository, rank(Hu) = Nu*Lw held while closed-loop performance ranged over two orders of magnitude, and cond(Hu) repeatedly ranked configurations backwards. The usable early warning is the magnitude of g: when median ‖g‖₁ climbs past a few hundred the configuration is heading for divergence whatever rank() reports.

The QP is condensed. Of the four equality blocks in the standard DeePC formulation, three are definitions — they pin u, y and σ_y to expressions in g, so substituting them out is exact, not an approximation. Only Up g = u_ini is a real restriction and stays a hard equality. That takes the problem from 287 variables and 208 equality rows to 91 and 12.

g is not sparse. Effective support is ~61 of 91 columns, and the prediction is built by large cancellation between opposing terms (weights of ±20 with ‖g‖₁ ≈ 500) rather than by selecting a few similar recorded trajectories. The regularizer on g here is an ℓ2 ridge, so that is expected — it has no mechanism to zero anything.

Solve time is the open problem. 156 ms median for rnd, 315 ms for lqr. Future work is the use OSQP to get faster solve times.