IMU Integration Error Simulator

Numerical Integration Method Comparison for IMU-Based Dead-Reckoning under Real Sensor Noise

A physics-based quadrotor simulator that demonstrates how integration method selection (Euler, Trapezoidal, Simpson) and IMU sensor noise interact to determine dead-reckoning position accuracy.

---

Overview

Drone flight controllers estimate position by double-integrating accelerometer readings — a technique called dead-reckoning. The accuracy depends on two factors: the numerical integration method chosen and the quality of the IMU sensor.

This project answers the question: does choosing a higher-order integration method (e.g., Simpson's rule over Euler) actually improve dead-reckoning accuracy in a real drone system?

The short answer is no — at typical IMU sampling rates (≥200 Hz), sensor noise dominates discretization error entirely. The noise floor is invariant to Δt, meaning switching from Euler to Simpson provides no measurable benefit under real MPU6050 noise conditions.

---

Reports

Full technical reports with derivations and results are available in both languages:

• English Report
• 한국어 보고서

---

Repository Structure

IMU-Simulator/
├── src/                  C simulation core
│   ├── main.c            Waypoint mission loop
│   ├── physics.c/h       RK4 quadrotor dynamics (12-state)
│   ├── controller.c/h    Cascaded PID controller
│   ├── imu_sim.c/h       MPU6050 noise model
│   ├── integrator.c/h    Euler / Trapezoidal / Simpson
│   └── logger.c/h        30-column CSV writer
├── scripts/
│   ├── visualize.py      Real-time 3D flight animation
│   └── analysis.py       Monte Carlo DT sweep + conclusions
├── docs/
│   ├── Report_EN.pdf
│   └── Report_KR.pdf
├── Makefile
└── run.sh                One-shot: build → simulate → visualize → analyze

---

How It Works

① Quadrotor flies a 10-waypoint mission under RK4 physics + PID control
② At each timestep, true acceleration is corrupted by MPU6050 noise model
③ Three integrators reconstruct position from the noisy acceleration
④ Logger records true vs. estimated positions for all three methods
⑤ visualize.py replays the flight in real-time 3D animation
⑥ analysis.py sweeps Δt from 1 ms to 100 ms across 80 Monte Carlo trials
⑦ Log-log convergence plot reveals the noise floor — all slopes ≈ 0

---

How to Run

Requirements: gcc, python3, numpy, matplotlib, pandas

bash run.sh

This builds the simulator, runs the 90-second waypoint mission, opens the 3D animation, then runs the DT sweep analysis and saves analysis.png.

Manual steps

make              # build
./imu_sim         # run simulation → result.csv
python3 scripts/visualize.py   # 3D animation
python3 scripts/analysis.py    # convergence analysis

---

Simulation

Quadrotor Dynamics

12-state Newton-Euler equations integrated with RK4 (Δt = 5 ms):

State	Symbol	Description
Position	x, y, z	Inertial frame (m)
Velocity	vx, vy, vz	Inertial frame (m/s)
Attitude	φ, θ, ψ	ZYX Euler angles (rad)
Angular rate	p, q, r	Body frame (rad/s)

Parameters

Parameter	Value
Mass	0.5 kg
Ixx = Iyy	5×10⁻³ kg·m²
Izz	9×10⁻³ kg·m²
Max thrust	9.81 N
Max torque	0.5 N·m

PID Controller

Cascaded architecture — outer position loop feeds the inner attitude loop:

Position error → [Kp_xy, Kd_xy, Ki_xy] → desired roll/pitch
              → [Kp_z,  Kd_z,  Ki_z ] → thrust
Attitude error → [Kp_att, Kd_att    ] → roll/pitch torques
               → [Kp_yaw, Kd_yaw    ] → yaw torque

Yaw reference is computed as atan2(Δy, Δx) toward each waypoint.

---

IMU Noise Model

MPU6050 datasheet values used for all noise sources:

Source	Model	Value
White noise	Gaussian, σ = NSD / √Δt	NSD = 400 μg/√Hz
Bias drift	Gauss-Markov, τ = 300 s	σ_b = 5 mg
Motor vibration	Gaussian	1.5 mg RMS
ADC quantization	Uniform rounding	±2g / 16-bit
Cross-axis sensitivity	Proportional coupling	2%

Gauss-Markov bias update at each step:

b_{k+1} = exp(−Δt/τ) · b_k  +  σ_b · √(1 − exp(−2Δt/τ)) · N(0,1)

---

Integration Methods

Position is reconstructed by double-integrating the noisy accelerometer output:

Method	Error order	FLOPs/step
Euler	O(Δt¹)	4
Trapezoidal	O(Δt²)	8
Simpson	O(Δt⁴)	10

Key Finding — Noise Floor Invariance

Position variance from white noise double-integration:

$$\sigma^2_{\text{pos}} = \frac{\text{NSD}^2 \cdot T^3}{3}$$

This is independent of Δt. At Δt = 5 ms, all three methods reach the same RMSE (~0.3 m over 10 s). The Trapezoidal method gives the best accuracy-per-FLOP under real sensor noise; Simpson's O(Δt⁴) advantage only appears in noiseless simulation.

---

Output CSV Format

result.csv — 30 columns logged at each timestep:

t, waypoint,
px, py, pz,          vx, vy, vz,
roll, pitch, yaw,    p, q, r,
ax_true, ay_true, az_true,
ax_noisy, ay_noisy, az_noisy,
px_euler, py_euler, pz_euler,
px_trap,  py_trap,  pz_trap,
px_simp,  py_simp,  pz_simp,
thrust

---

License

Licensed under the Apache License 2.0.