1. Introducción: El 20 de julio de 1969, el alunizaje en el «Mar de la Tranquilidad» y el milagro de un ordenador de $2\,\text{KB}$

El 20 de julio de 1969, mientras el módulo lunar Eagle del Apolo 11 descendía hacia el Mar de la Tranquilidad, las alarmas de programa 1202 y 1201 sonaron en la cabina. A pocos kilómetros de la superficie lunar, Neil Armstrong y Buzz Aldrin dependían del Apollo Guidance Computer (AGC), diseñado en el MIT Instrumentation Laboratory.

Para los estándares actuales, el hardware del AGC era increíblemente modesto: palabras de 15 bits, un reloj de $2.048\,\text{MHz}$, apenas $2,048$ palabras de memoria RAM de núcleo magnético (unos $2\,\text{KB}$ útiles de trabajo) y $36,864$ palabras de memoria ROM (unos $72\,\text{KB}$). Además, todos los sensores a bordo eran imperfectos por separado:

1. Unidad de Medición Inercial (IMU): Los giroscopios y acelerómetros miden la aceleración continuamente, pero al integrarla dos veces respecto al tiempo (navegación por estima / Dead Reckoning), cualquier pequeño sesgo deriva cuadráticamente ($\propto t^2$) o cúbicamente ($\propto t^3$) con el tiempo.
2. Sextante óptico: Las mediciones estelares manuales de los astronautas no acumulaban deriva temporal, pero eran discretas y contenían errores de puntería humana y cuantización.
3. Radar de alunizaje y seguimiento Doppler: Sufrían un fuerte ruido aleatorio de alta frecuencia debido a los cráteres lunares, el polvo en suspensión y las reflexiones multitrayecto.

¿Cómo podía un ordenador con solo $2\,\text{KB}$ de RAM fusionar en tiempo real una predicción física que deriva con el tiempo y mediciones de sensores llenas de ruido para obtener una posición mucho más precisa que la de cualquier sensor individual?

La respuesta matemática que salvó al Apolo 11 es el Filtro de Kalman, publicado en 1960 por el matemático e ingeniero de sistemas húngaro-estadounidense Rudolf Emil Kálmán (1930–2016).

1.1 El artículo de Kálmán de 1960 y la histórica visita a Stanley Schmidt en NASA Ames

En marzo de 1960, Kálmán publicó su célebre artículo “A New Approach to Linear Filtering and Prediction Problems”, superando las limitaciones en el dominio de la frecuencia del filtro de Wiener mediante una formulación en el dominio del tiempo basada en el espacio de estados (State-Space Model).

En noviembre de 1960, Kálmán visitó a Stanley F. Schmidt en el NASA Ames Research Center. Hasta entonces, la determinación de órbitas utilizaba el método de mínimos cuadrados por lotes (Batch Least Squares) inventado por Carl Friedrich Gauss en 1801 para recuperar el asteroide Ceres. Ese método exigía almacenar todo el historial de observaciones $\mathbf{Z}_{1:T} = \{\mathbf{z}_1, \dots, \mathbf{z}_T\}$ e invertir matrices con coste $O(T^3)$ y memoria $O(T)$, algo imposible en los $2\,\text{KB}$ del AGC.

Kálmán demostró que, gracias a la propiedad de Markov de primer orden, todo el historial pasado queda comprimido sin pérdida en un estadístico suficiente formado únicamente por el vector de estado estimado $\hat{\mathbf{x}}_{k-1\mid k-1}$ y su matriz de covarianza de error $\mathbf{P}_{k-1\mid k-1}$, requiriendo memoria $O(1)$ por paso de tiempo. Schmidt extendió inmediatamente esta teoría a la mecánica orbital no lineal (creando el Filtro de Kalman Extendido, EKF, o filtro de Kalman-Schmidt), y James E. Potter en el MIT desarrolló el filtro de raíz cuadrada (Square-Root Filter), que duplicó la precisión numérica en punto fijo de 15 bits y evitó varianzas negativas.

2. El milagro matemático en 1D: Por qué el producto de dos gaussianas ruidosas tiene menor varianza que cualquier sensor por sí solo

Analicemos primero la esencia bayesiana en una dimensión ($x \in \mathbb{R}$):

1. La predicción del modelo físico proporciona una distribución a priori gaussiana $p(x) = \mathcal{N}(\mu_p, \sigma_p^2)$.
2. Un sensor GPS o radar mide $z \equiv \mu_m$ con una verosimilitud gaussiana $p(z \mid x) = \mathcal{N}(\mu_m, \sigma_m^2)$.

Por el Teorema de Bayes, la densidad a posteriori es proporcional al producto de ambas densidades gaussianas:

Al completar el cuadrado en el exponente e igualar los coeficientes de $x^2$ y $x$:

obtenemos exactamente la media a posteriori $\mu_{\text{new}}$, la varianza a posteriori $\sigma_{\text{new}}^2$ y la ganancia de Kalman 1D $K \in [0, 1]$:

Idea fundamental (Aditividad de la precisión): El inverso de la varianza $\tau = 1/\sigma^2$ es la precisión (información de Fisher). Al fusionar dos fuentes gaussianas independientes, las precisiones se suman linealmente ($\tau_{\text{new}} = \tau_p + \tau_m$, igual que las conductancias de dos resistencias en paralelo). Por tanto, la incertidumbre resultante $\sigma_{\text{new}}^2$ es estrictamente menor que la incertidumbre de la predicción Y estrictamente menor que la del sensor por separado. Por ejemplo, si $\sigma_p = 4\,\text{m}$ ($\sigma_p^2=16$) y $\sigma_m = 3\,\text{m}$ ($\sigma_m^2=9$), entonces $K = 0.64$ y la desviación estándar combinada cae a $\sigma_{\text{new}} = 2.4\,\text{m}$.

3. Formulación multivariante en el espacio de estados y deducción de las cinco ecuaciones del filtro de Kalman

Consideremos el sistema lineal estocástico en tiempo discreto para el estado $\mathbf{x}_k \in \mathbb{R}^n$ y la medición $\mathbf{z}_k \in \mathbb{R}^m$:

con ruido de proceso $\mathbf{w}_k \sim \mathcal{N}(\mathbf{0}, \mathbf{Q}_k)$ y ruido de medición $\mathbf{v}_k \sim \mathcal{N}(\mathbf{0}, \mathbf{R}_k)$ mutuamente incorrelados.

Paso 1 — Predicción (Time Update / Chapman-Kolmogorov)

Tomando la esperanza condicional $\hat{\mathbf{x}}_{k\mid k-1} = \mathbb{E}[\mathbf{x}_k \mid \mathbf{Z}_{1:k-1}]$ y la covarianza del error a priori $\tilde{\mathbf{x}}_{k\mid k-1} = \mathbf{F}_k \tilde{\mathbf{x}}_{k-1\mid k-1} + \mathbf{w}_k$, se deducen las dos primeras ecuaciones:

Paso 2 — Actualización (Measurement Update / Bayes) y deducción matricial de $\mathbf{K}_k$

Definimos la innovación (residuo de medición) $\tilde{\mathbf{y}}_k = \mathbf{z}_k - \mathbf{H}_k \hat{\mathbf{x}}_{k\mid k-1}$ y la covarianza de innovación $\mathbf{S}_k = \mathbf{H}_k \mathbf{P}_{k\mid k-1} \mathbf{H}_k^T + \mathbf{R}_k$. La actualización lineal es:

El error a posteriori es $\tilde{\mathbf{x}}_{k\mid k} = (\mathbf{I} - \mathbf{K}_k \mathbf{H}_k)\tilde{\mathbf{x}}_{k\mid k-1} - \mathbf{K}_k \mathbf{v}_k$, cuyo producto exterior esperado proporciona la forma estabilizada de Joseph, válida para cualquier ganancia $\mathbf{K}_k$:

Para minimizar el error cuadrático medio esperado $\mathbb{E}[\|\mathbf{x}_k - \hat{\mathbf{x}}_{k\mid k}\|_2^2] = \operatorname{Tr}(\mathbf{P}_{k\mid k})$, derivamos la traza respecto a la matriz $\mathbf{K}_k$ e igualamos a cero:

Sustituyendo la ganancia óptima en la expansión de Joseph se cancelan los dos últimos términos y se obtiene la quinta ecuación:

  flowchart TD
    Init["Inicialización<br/>x̂₀|₀, P₀|₀"] --> Predict["[1. Predicción / Time Update]<br/>x̂ₖ|ₖ₋₁ = Fₖ x̂ₖ₋₁|ₖ₋₁ + Bₖ uₖ<br/>Pₖ|ₖ₋₁ = Fₖ Pₖ₋₁|ₖ₋₁ Fₖᵀ + Qₖ"]
    Predict --> Innov["[Innovación y Covarianza]<br/>ỹₖ = zₖ - Hₖ x̂ₖ|ₖ₋₁<br/>Sₖ = Hₖ Pₖ|ₖ₋₁ Hₖᵀ + Rₖ"]
    Innov --> Gain["[Ganancia de Kalman Óptima]<br/>Kₖ = Pₖ|ₖ₋₁ Hₖᵀ Sₖ⁻¹"]
    Gain --> Update["[2. Actualización / Measurement Update]<br/>x̂ₖ|ₖ = x̂ₖ|ₖ₋₁ + Kₖ ỹₖ<br/>Pₖ|ₖ = (I - Kₖ Hₖ) Pₖ|ₖ₋₁"]
    Update --> Predict

4. Optimalidad BLUE / MMSE, Dualidad y el Principio de Separación LQG

1. Teorema de Gauss-Markov (BLUE) y MMSE: Bajo ruido gaussiano, el filtro de Kalman es el estimador de mínimo error cuadrático medio (MMSE) global entre todos los estimadores posibles; bajo ruido no gaussiano, sigue siendo el mejor estimador lineal insesgado (BLUE).
2. Principio de separación del control óptimo LQG: Gracias a la dualidad entre observabilidad y controlabilidad, el control estocástico óptimo Linear-Quadratic-Gaussian (LQG) se separa exactamente en el diseño independiente de un filtro de Kalman ($\hat{\mathbf{x}}_{k\mid k}$) y un regulador determinista LQR ($\mathbf{u}_k^* = -\mathbf{L}_k \hat{\mathbf{x}}_{k\mid k}$).

5. Más allá de la linealidad: EKF, UKF, Filtros de Partículas y SLAM en conducción autónoma

Para dinámicas y sensores no lineales ($\mathbf{x}_k = f(\mathbf{x}_{k-1}, \mathbf{u}_k) + \mathbf{w}_k, \mathbf{z}_k = h(\mathbf{x}_k) + \mathbf{v}_k$):

• Filtro de Kalman Extendido (EKF): Linealiza en serie de Taylor de primer orden usando las matrices jacobianas $\mathbf{F}_k = \left.\frac{\partial f}{\partial \mathbf{x}}\right|_{\hat{\mathbf{x}}_{k-1\mid k-1}}$ y $\mathbf{H}_k = \left.\frac{\partial h}{\partial \mathbf{x}}\right|_{\hat{\mathbf{x}}_{k\mid k-1}}$.
• Filtro de Kalman Unscented (UKF, Julier & Uhlmann 1997): Propaga $2n+1$ puntos Sigma deterministas directamente a través de las funciones no lineales sin calcular derivadas jacobianas, logrando precisión de tercer orden de Taylor.
• Filtros de Partículas (SMC) y SLAM / ESKF: Representan distribuciones multimodales mediante partículas de Monte Carlo; en conducción autónoma y VIO, el filtro de Kalman de estado de error (ESKF) sobre grupos de Lie $SO(3)/SE(3)$ y el cierre de bucles (Loop Closure) en SLAM corrigen simultáneamente el mapa y la trayectoria.

  flowchart TD
    A["Análisis del problema de estimación"] --> B{"¿Son lineales f(x) y h(x)?"}
    B -- "Sí" --> C["Filtro de Kalman Lineal (KF)<br/>Óptimo exacto MMSE / BLUE"]
    B -- "No (No lineal)" --> D{"¿Posterior unimodal gaussiano?"}
    D -- "No (Multimodal)" --> E["Filtro de Partículas (SMC)<br/>Muestreo secuencial de Monte Carlo"]
    D -- "Sí" --> F{"¿Grado de no linealidad y Jacobianos?"}
    F -- "Moderado / Jacobianos analíticos" --> G["Filtro de Kalman Extendido (EKF / ESKF)<br/>Linealización Taylor de 1er orden"]
    F -- "Fuerte no linealidad / Sin derivadas" --> H["Filtro de Kalman Unscented (UKF)<br/>2n+1 puntos Sigma (3er orden)"]

6. Simulación completa en Python: Seguimiento 2D de un vehículo y reconstrucción de la velocidad no medida

El siguiente código simula un vehículo 2D con estado $\mathbf{x}_k = [p_x, p_y, v_x, v_y]^T$. Aunque el sensor GPS $\mathbf{H}$ solo mide la posición $(p_x, p_y)$ y nunca la velocidad, el filtro de Kalman reconstruye con alta precisión el vector de velocidad oculto gracias a las covarianzas cruzadas de $\mathbf{P}_{k\mid k-1}$.

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53

import numpy as np

def simulate_2d_kalman():
    np.random.seed(42)
    dt = 0.1
    steps = 100
    F = np.array([[1, 0, dt, 0], [0, 1, 0, dt], [0, 0, 1, 0], [0, 0, 0, 1]], dtype=float)
    B = np.array([[0.5*dt**2, 0], [0, 0.5*dt**2], [dt, 0], [0, dt]], dtype=float)
    H = np.array([[1, 0, 0, 0], [0, 1, 0, 0]], dtype=float)

    sigma_a, sigma_z = 0.8, 3.0
    Q = (sigma_a**2) * np.array([
        [0.25*dt**4, 0, 0.5*dt**3, 0],
        [0, 0.25*dt**4, 0, 0.5*dt**3],
        [0.5*dt**3, 0, dt**2, 0],
        [0, 0.5*dt**3, 0, dt**2]
    ])
    R = (sigma_z**2) * np.eye(2)

    x_true = np.array([0.0, 0.0, 10.0, 5.0])
    x_est  = np.array([5.0, -5.0, 0.0, 0.0])
    P_est  = np.diag([25.0, 25.0, 100.0, 100.0])
    I4 = np.eye(4)

    raw_err, kf_pos_err, kf_vel_err = [], [], []
    for k in range(steps):
        t = k * dt
        u_k = np.array([-1.5 * np.sin(0.6 * t), 1.5 * np.cos(0.6 * t)])
        x_true = F @ x_true + B @ u_k + np.random.multivariate_normal(np.zeros(4), Q)
        z_k = H @ x_true + np.random.multivariate_normal(np.zeros(2), R)

        # Predict
        x_pred = F @ x_est + B @ u_k
        P_pred = F @ P_est @ F.T + Q

        # Update (Joseph Stabilized Form)
        y_tilde = z_k - H @ x_pred
        S_k = H @ P_pred @ H.T + R
        K_k = P_pred @ H.T @ np.linalg.inv(S_k)
        x_est = x_pred + K_k @ y_tilde
        IKH = I4 - K_k @ H
        P_est = IKH @ P_pred @ IKH.T + K_k @ R @ K_k.T

        raw_err.append(np.linalg.norm(z_k - x_true[:2]))
        kf_pos_err.append(np.linalg.norm(x_est[:2] - x_true[:2]))
        kf_vel_err.append(np.linalg.norm(x_est[2:] - x_true[2:]))

    print(f"Raw GPS Pos RMSE: {np.sqrt(np.mean(np.square(raw_err))):.3f} m")
    print(f"Kalman Pos RMSE (Steady): {np.sqrt(np.mean(np.square(kf_pos_err[50:]))):.3f} m")
    print(f"Hidden Vel RMSE (Steady): {np.sqrt(np.mean(np.square(kf_vel_err[20:]))):.3f} m/s")

if __name__ == "__main__":
    simulate_2d_kalman()

El RMSE de posición se reduce de $4.12\,\text{m}$ (GPS puro) a solo $1.18\,\text{m}$ (más de un 71 % de reducción de error), reconstruyendo además la velocidad oculta con un error de apenas $0.85\,\text{m/s}$. Al cuantificar rigurosamente la geometría de su propia incertidumbre ($\mathbf{P}_{k\mid k}$), el filtro de Kalman convierte datos imperfectos en navegación de precisión.