Revelando el poder del filtro de Kalman

Descubriendo el poder del filtro de Kalman

Como científicos de datos, a veces nos enfrentamos a situaciones en las que necesitamos modelar una tendencia para predecir valores futuros. Aunque hay una tentación de centrarse en algoritmos estadísticos o basados en aprendizaje automático, estoy aquí para presentar una opción diferente: el Filtro de Kalman (KF).

En la década de 1960, Rudolf E. Kalman revolucionó cómo se pueden modelar sistemas complejos con el KF. Desde guiar aeronaves o naves espaciales hacia su destino hasta permitir que tu teléfono inteligente encuentre su lugar en este mundo, este algoritmo combina datos y matemáticas para proporcionar estimaciones de futuros estados con una precisión increíble.

En este blog profundizaremos en cómo funciona el Filtro de Kalman, mostrando ejemplos en Python que enfatizarán el verdadero poder de esta técnica. Comenzando con un ejemplo simple en 2D, veremos cómo podemos modificar el código para adaptarlo a espacios más avanzados en 4D y terminaremos cubriendo el Filtro de Kalman Extendido (el sofisticado predecesor). Únete a mí en este viaje mientras nos adentramos en el mundo de los algoritmos y filtros predictivos.

Los fundamentos de un filtro de Kalman

El KF proporciona una estimación del estado de un sistema construyendo y actualizando continuamente un conjunto de matrices de covarianza (que representan la distribución estadística de ruido y estados pasados) recopiladas a partir de observaciones y otras medidas en el tiempo. A diferencia de otros algoritmos listos para usar, es posible expandir y refinar directamente la solución definiendo las relaciones matemáticas entre el sistema y las fuentes externas. Aunque esto puede sonar bastante complejo e intrincado, el proceso se puede resumir en dos pasos: predecir y actualizar. Estas fases trabajan juntas para corregir y refinar de manera iterativa las estimaciones de estado del sistema.

Paso de predicción:

Esta fase se trata de pronosticar el próximo estado futuro del sistema en función de las estimaciones a posteriori conocidas del modelo y el paso en el tiempo de Δk. Matemáticamente, representamos las estimaciones del espacio de estados como:

donde F, nuestra matriz de transición de estados, modela cómo se desarrollan los estados de un paso a otro, independientemente de la entrada de control y del ruido del proceso. Nuestra matriz B modela la influencia de la entrada de control, uₖ, en el estado.

Junto con nuestras estimaciones del próximo estado, el algoritmo también calcula la incertidumbre de la estimación representada por la matriz de covarianza P:

La matriz de covarianza predictiva del estado representa la confianza y precisión de nuestras predicciones, influenciada por Q, la matriz de covarianza del ruido del proceso del propio sistema. Aplicamos esta matriz a ecuaciones subsiguientes en el paso de actualización para corregir la información que el Filtro de Kalman tiene sobre el sistema, mejorando así las estimaciones de los futuros estados.

Paso de actualización:

En el paso de actualización, el algoritmo realiza actualizaciones en la Ganancia de Kalman, las estimaciones de estado y la matriz de covarianza. La Ganancia de Kalman determina cuánta influencia debería tener una nueva medición en las estimaciones de estado. El cálculo incluye la matriz del modelo de observación, H, que relaciona el estado con la medición que esperamos recibir, y R, la matriz de covarianza del ruido de la medición de los errores:

En esencia, K intenta equilibrar la incertidumbre en las predicciones con la presente en las mediciones. Como se mencionó anteriormente, la Ganancia de Kalman se aplica para corregir las estimaciones de estado y la covarianza, como se presentan en las siguientes ecuaciones respectivamente:

donde el cálculo entre corchetes para la estimación de estado es el residual, la diferencia entre la medición real y la que el modelo predijo.

La verdadera belleza de cómo funciona un Filtro de Kalman radica en su naturaleza recursiva, actualizando continuamente tanto el estado como la covarianza a medida que se recibe nueva información. Esto permite que el modelo se refine de manera estadísticamente óptima con el tiempo, lo cual es un enfoque particularmente poderoso para modelar sistemas que reciben una corriente de observaciones ruidosas.

El filtro de Kalman en acción

Es posible sentirse abrumado por las ecuaciones que sustentan el Filtro de Kalman, y para apreciar completamente cómo funciona solo con las matemáticas, sería necesario entender el espacio de estado (fuera del alcance de este blog), pero intentaré darle vida con algunos ejemplos en Python. En su forma más simple, podemos definir un objeto Filtro de Kalman como:

import numpy as npclass KalmanFilter:    """    Una implementación del clásico Filtro de Kalman para sistemas dinámicos lineales.    El Filtro de Kalman es un algoritmo óptimo de procesamiento de datos recursivo que    tiene como objetivo estimar el estado de un sistema a partir de observaciones ruidosas.    Atributos:        F (np.ndarray): La matriz de transición de estado.        B (np.ndarray): La matriz de entrada de control.        H (np.ndarray): La matriz de observación.        u (np.ndarray): la entrada de control.        Q (np.ndarray): La matriz de covarianza de ruido de proceso.        R (np.ndarray): La matriz de covarianza de ruido de medición.        x (np.ndarray): La estimación de estado media del paso anterior (k-1).        P (np.ndarray): La covarianza de estado del paso anterior (k-1).    """    def __init__(self, F, B, u, H, Q, R, x0, P0):        """        Inicializa el Filtro de Kalman con las matrices necesarias y el estado inicial.        Parámetros:            F (np.ndarray): La matriz de transición de estado.            B (np.ndarray): La matriz de entrada de control.            H (np.ndarray): La matriz de observación.            u (np.ndarray): la entrada de control.            Q (np.ndarray): La matriz de covarianza de ruido de proceso.            R (np.ndarray): La matriz de covarianza de ruido de medición.            x0 (np.ndarray): La estimación de estado inicial.            P0 (np.ndarray): La covarianza de estado inicial.        """        self.F = F  # Matriz de transición de estado        self.B = B  # Matriz de entrada de control        self.u = u  # Vector de control        self.H = H  # Matriz de observación        self.Q = Q  # Covarianza de ruido de proceso        self.R = R  # Covarianza de ruido de medición        self.x = x0  # Estimación de estado inicial        self.P = P0  # Covarianza de estimación inicial    def predict(self):        """        Predice el estado y la covarianza del estado para el siguiente paso en el tiempo.        """        self.x = self.F @ self.x + self.B @ self.u        self.P = self.F @ self.P @ self.F.T + self.Q        return self.x    def update(self, z):        """        Actualiza la estimación de estado con la última medición.        Parámetros:            z (np.ndarray): La medición en el paso actual.        """        y = z - self.H @ self.x        S = self.H @ self.P @ self.H.T + self.R        K = self.P @ self.H.T @ np.linalg.inv(S)        self.x = self.x + K @ y        I = np.eye(self.P.shape[0])        self.P = (I - K @ self.H) @ self.P                return self.xDesafíos con sistemas no lineales

Utilizaremos las funciones predict() y update() para iterar a través de los pasos descritos anteriormente. El diseño del filtro anterior funcionará para cualquier serie de tiempo, y para mostrar cómo nuestras estimaciones se comparan con los valores reales, construyamos un ejemplo simple:

import numpy as npimport matplotlib.pyplot as plt# Establecer la semilla aleatoria para reproducibilidadnp.random.seed(42)# Simular la posición real del objeto con el tiempo.true_velocity = 0.5  # unidades por paso de tiemponum_steps = 50time_steps = np.linspace(0, num_steps, num_steps)true_positions = true_velocity * time_steps# Simular las mediciones con ruidomeasurement_noise = 10  # aumentar este valor para hacer que las mediciones sean más ruidosasnoisy_measurements = true_positions + np.random.normal(0, measurement_noise, num_steps)# Graficar las posiciones reales y las mediciones ruidosasplt.figure(figsize=(10, 6))plt.plot(time_steps, true_positions, label='Posición Real', color='green')plt.scatter(time_steps, noisy_measurements, label='Mediciones Ruidosas', color='red', marker='o')plt.xlabel('Paso de Tiempo')plt.ylabel('Posición')plt.title('Posición Real y Mediciones Ruidosas a lo largo del tiempo')plt.legend()plt.show()

En realidad, la ‘Posición Real’ sería desconocida, pero la graficaremos aquí como referencia. Las ‘Mediciones Ruidosas’ son los puntos de observación que se alimentan en nuestro Filtro de Kalman. Realizaremos una instanciación muy básica de las matrices, y en cierta medida no importa ya que el modelo de Kalman convergerá rápidamente mediante la aplicación de la Ganancia de Kalman, pero podría ser razonable en ciertas circunstancias realizar un arranque en caliente del modelo.

# Inicialización del Filtro de Kalman
F = np.array([[1, 1], [0, 1]])   # Matriz de transición de estados
B = np.array([[0], [0]])          # Sin entrada de control
u = np.array([[0]])               # Sin entrada de control
H = np.array([[1, 0]])            # Función de medición
Q = np.array([[1, 0], [0, 3]])    # Matriz de covarianza de ruido del proceso
R = np.array([[measurement_noise**2]]) # Matriz de covarianza de ruido de medición
x0 = np.array([[0], [0]])         # Estimación inicial de estado
P0 = np.array([[1, 0], [0, 1]])   # Covarianza de estimación inicial
kf = KalmanFilter(F, B, u, H, Q, R, x0, P0)

# Espacio para posiciones y velocidades estimadas
estimated_positions = np.zeros(num_steps)
estimated_velocities = np.zeros(num_steps)

# Bucle del Filtro de Kalman
for t in range(num_steps):
    # Predecir
    kf.predict()
    
    # Actualizar
    measurement = np.array([[noisy_measurements[t]]])
    kf.update(measurement)
    
    # Almacenar la posición y velocidad filtrada
    estimated_positions[t] = kf.x[0]
    estimated_velocities[t] = kf.x[1]

# Graficar las posiciones verdaderas, las mediciones ruidosas y las estimaciones del Filtro de Kalman
plt.figure(figsize=(10, 6))
plt.plot(time_steps, true_positions, label='Posición Verdadera', color='green')
plt.scatter(time_steps, noisy_measurements, label='Mediciones Ruidosas', color='red', marker='o')
plt.plot(time_steps, estimated_positions, label='Estimación del Filtro de Kalman', color='blue')
plt.xlabel('Paso de Tiempo')
plt.ylabel('Posición')
plt.title('Estimación del Filtro de Kalman a lo largo del Tiempo')
plt.legend()
plt.show()

Incluso con este diseño muy simple de una solución, el modelo hace un trabajo aceptable al atravesar el ruido para encontrar la posición verdadera. Esto podría funcionar bien para aplicaciones simples, pero a menudo las tendencias son más sutiles y se ven afectadas por eventos externos. Para manejar esto, generalmente necesitamos modificar tanto la representación del espacio de estados como la forma en que calculamos las estimaciones y corregimos la matriz de covarianza cuando llega nueva información, vamos a explorar esto más con otro ejemplo.

Rastreando un objeto en movimiento en 4D

Supongamos que queremos rastrear el movimiento de un objeto en el espacio y el tiempo, y para hacer este ejemplo aún más realista, simularemos alguna fuerza que actúa sobre él y resulta en una rotación angular. Para mostrar la adaptabilidad de este algoritmo a representaciones de espacios de estados de mayor dimensión, asumiremos una fuerza lineal, aunque en realidad este no será el caso (exploraremos un ejemplo más realista después de esto). El siguiente código muestra un ejemplo de cómo modificaríamos nuestro Filtro de Kalman para este escenario en particular.

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

class KalmanFilter:
"""Una implementación del clásico Filtro de Kalman para sistemas dinámicos lineales.
El Filtro de Kalman es un algoritmo óptimo de procesamiento de datos recursivos que
tiene como objetivo estimar el estado de un sistema a partir de observaciones ruidosas.

Atributos:
F (np.ndarray): Matriz de transición de estados.
B (np.ndarray): Matriz de entrada de control.
H (np.ndarray): Matriz de observación.
u (np.ndarray): Entrada de control.
Q (np.ndarray): Matriz de covarianza de ruido del proceso.
R (np.ndarray): Matriz de coviarianza de ruido de medición.
x (np.ndarray): Estimación media del estado del paso anterior (k-1).
P (np.ndarray): Covarianza de estimación del paso anterior (k-1).
"""

def __init__(self, F=None, B=None, u=None, H=None, Q=None, R=None, x0=None, P0=None):
"""Inicializa el Filtro de Kalman con las matrices necesarias y el estado inicial.

Parámetros:
F (np.ndarray): Matriz de transición de estados.
B (np.ndarray): Matriz de entrada de control.
H (np.ndarray): Matriz de observación.
u (np.ndarray): Entrada de control.
Q (np.ndarray): Matriz de covarianza de ruido del proceso.
R (np.ndarray): Matriz de coviarianza de ruido de medición.
x0 (np.ndarray): Estimación inicial del estado.
P0 (np.ndarray): Matriz de covarianza de estimación inicial.
"""

self.F = F # Matriz de transición de estados
self.B = B # Matriz de entrada de control
self.u = u # Entrada de control
self.H = H # Matriz de observación
self.Q = Q # Matriz de covarianza de ruido del proceso
self.R = R # Matriz de covarianza de ruido de medición
self.x = x0 # Estimación inicial de estado
self.P = P0 # Covarianza de estimación inicial

def predict(self):
"""Predice el estado y la covarianza del estado para el próximo paso de tiempo."""

self.x = np.dot(self.F, self.x) + np.dot(self.B, self.u)
self.P = np.dot(np.dot(self.F, self.P), self.F.T) + self.Q

def update(self, z):
"""Actualiza la estimación del estado con la última medición.

Parámetros:
z (np.ndarray): La medición en el paso actual.
"""

y = z - np.dot(self.H, self.x)
S = np.dot(self.H, np.dot(self.P, self.H.T)) + self.R
K = np.dot(np.dot(self.P, self.H.T), np.linalg.inv(S))
self.x = self.x + np.dot(K, y)
self.P = self.P - np.dot(np.dot(K, self.H), self.P)

# Parámetros para la simulación
true_angular_velocity = 0.1 # radianes por paso de tiempo
radius = 20
num_steps = 100
dt = 1 # paso de tiempo

# Crear pasos de tiempo
time_steps = np.arange(0, num_steps*dt, dt)

# Estado verdadero del suelo
true_x_positions = radius * np.cos(true_angular_velocity * time_steps)
true_y_positions = radius * np.sin(true_angular_velocity * time_steps)
true_z_positions = 0.5 * time_steps # velocidad constante en z

# Crear mediciones ruidosas
measurement_noise = 1.0
noisy_x_measurements = true_x_positions + np.random.normal(0, measurement_noise, num_steps)
noisy_y_measurements = true_y_positions + np.random.normal(0, measurement_noise, num_steps)
noisy_z_measurements = true_z_positions + np.random.normal(0, measurement_noise, num_steps)

# Inicialización del Filtro de Kalman
F = np.array([[1, 0, 0, -radius*dt*np.sin(true_angular_velocity*dt)],
[0, 1, 0, radius*dt*np.cos(true_angular_velocity*dt)],
[0, 0, 1, 0],
[0, 0, 0, 1]])
B = np.zeros((4, 1))
u = np.zeros((1, 1))
H = np.array([[1, 0, 0, 0],
[0, 1, 0, 0],
[0, 0, 1, 0]])
Q = np.eye(4) * 0.1 # Ruido de proceso pequeño
R = measurement_noise**2 * np.eye(3) # Ruido de medición
x0 = np.array([[0], [radius], [0], [true_angular_velocity]])
P0 = np.eye(4) * 1.0
kf = KalmanFilter(F, B, u, H, Q, R, x0, P0)

#

Algunos puntos interesantes para señalar aquí, en el gráfico anterior podemos ver cómo el modelo se corrige rápidamente al estado verdadero estimado a medida que comenzamos a iterar sobre las observaciones. El modelo también funciona bastante bien al identificar el estado verdadero del sistema, con estimaciones que se cruzan con los estados verdaderos ('trayectoria verdadera'). Este diseño podría ser apropiado para algunas aplicaciones del mundo real, pero no para aquellas en las que la fuerza que actúa sobre el sistema es no lineal. En su lugar, debemos considerar una aplicación diferente del Filtro de Kalman: el Filtro de Kalman Extendido, un precursor de lo que hemos explorado hasta ahora que lineariza la no linealidad de las observaciones entrantes.

El Filtro de Kalman Extendido

Cuando intentamos modelar un sistema en el que las observaciones o la dinámica del sistema son no lineales, debemos aplicar el Filtro de Kalman Extendido (EKF). Este algoritmo difiere del anterior al incorporar la matriz Jacobiana en la solución y realizar una expansión de la serie de Taylor para encontrar aproximaciones lineales de primer orden de las modelos de transición y observación del estado. Para expresar esta extensión matemáticamente, nuestros cálculos algorítmicos clave ahora son:

para la predicción del estado, donde f es nuestra función de transición de estado no lineal aplicada a la estimación del estado anterior y u es la entrada de control en el paso de tiempo anterior.

para la predicción de la covarianza del error, donde F es la Jacobiana de la función de transición de estado con respecto a P la covarianza del error anterior y Q la matriz de covarianza del ruido del proceso.

la observación de nuestra medición, z, en el paso de tiempo k, donde h es la función de observación no lineal aplicada a nuestra predicción de estado con cierto ruido de observación aditivo v.

nuestra actualización del cálculo de la ganancia de Kalman, con H la Jacobiana de la función de observación con respecto al estado y R la matriz de covarianza del ruido de medición.

el cálculo modificado para la estimación del estado que incorpora la ganancia de Kalman y la función de observación no lineal, y finalmente nuestra ecuación para actualizar la covarianza del error:

En el último ejemplo, esto utilizará la Jacobiana para linealizar el efecto no lineal de la rotación angular en nuestro objeto, modificando el código adecuadamente. Diseñar un EKF es más desafiante que el KF, ya que nuestra suposición de aproximaciones lineales de primer orden puede introducir inadvertidamente errores en nuestras estimaciones de estado. Además, el cálculo de la Jacobiana puede volverse complejo, computacionalmente costoso y difícil de definir para ciertas situaciones, lo que también puede conducir a errores. Sin embargo, si se diseña correctamente, el EKF a menudo superará la implementación del KF.

Basándonos en nuestro último ejemplo en Python, he presentado la implementación del EKF:

import numpy as npimport matplotlib.pyplot as pltfrom mpl_toolkits.mplot3d import Axes3Dclass ExtendedKalmanFilter:    """    Una implementación del Filtro de Kalman Extendido (EKF).    Este filtro es adecuado para sistemas con dinámica no lineal al linealizar    el modelo del sistema en cada paso de tiempo utilizando la matriz Jacobiana.        Atributos:        state_transition (callable): La función de transición de estado para el sistema.        jacobian_F (callable): Función para calcular la matriz Jacobiana de la función de transición de estado.        H (np.ndarray): La matriz de observación.        jacobian_H (callable): Función para calcular la matriz Jacobiana del modelo de observación.        Q (np.ndarray): La matriz de covarianza del ruido de proceso.        R (np.ndarray): La matriz de covarianza del ruido de medición.        x (np.ndarray): La estimación inicial del estado.        P (np.ndarray): La covarianza inicial de la estimación.    """    def __init__(self, state_transition, jacobian_F, observation_matrix, jacobian_H, Q, R, x, P):        """        Construye el Filtro de Kalman Extendido.        Parámetros:            state_transition (callable): La función de transición de estado.            jacobian_F (callable): Función para calcular la matriz Jacobiana de F.            observation_matrix (np.ndarray): Matriz de observación.            jacobian_H (callable): Función para calcular la matriz Jacobiana de H.            Q (np.ndarray): Covarianza del ruido de proceso.            R (np.ndarray): Covarianza del ruido de medición.            x (np.ndarray): Estimación inicial del estado.            P (np.ndarray): Covarianza inicial de la estimación.        """        self.state_transition = state_transition  # Función de transición de estado no lineal        self.jacobian_F = j

Un resumen sencillo

En este blog hemos explorado en profundidad cómo construir y aplicar un Filtro de Kalman (KF), así como también cómo implementar un Filtro de Kalman Extendido (EKF). Terminemos resumiendo los casos de uso, ventajas y desventajas de estos modelos.

KF: Este modelo es aplicable a sistemas lineales, donde podemos asumir que las transiciones de estado y las matrices de observación son funciones lineales del estado (con cierto ruido gaussiano). Podrías considerar aplicar este algoritmo para:

  • seguimiento de la posición y velocidad de un objeto que se mueve a una velocidad constante
  • aplicaciones de procesamiento de señales si el ruido es estocástico o se puede representar mediante modelos lineales
  • previsión económica si las relaciones subyacentes se pueden modelar de forma lineal

La principal ventaja del KF es que (una vez que sigas los cálculos matriciales) el algoritmo es bastante simple, menos intensivo computacionalmente que otros enfoques y puede proporcionar pronósticos y estimaciones muy precisos del estado real del sistema en el tiempo. La desventaja es la suposición de linealidad, que típicamente no se cumple en escenarios complejos del mundo real.

EKF: Podemos considerar esencialmente al EKF como el equivalente no lineal del KF, respaldado por el uso del Jacobiano. Considerarías el EKF si estás tratando con:

  • sistemas robóticos donde las medidas y la dinámica del sistema suelen ser no lineales
  • sistemas de seguimiento y navegación que a menudo involucran velocidades no constantes o cambios en la velocidad angular como el seguimiento de aeronaves o naves espaciales.
  • sistemas automotrices al implementar control de crucero o mantenimiento de carril que se encuentra en los autos "inteligentes" más modernos.

El EKF a menudo produce mejores estimaciones que el KF, especialmente para sistemas no lineales, pero puede volverse mucho más intensivo computacionalmente debido al cálculo de la matriz Jacobiana. Este enfoque también se basa en aproximaciones lineales de primer orden de la expansión de la serie de Taylor, lo cual puede no ser una suposición apropiada para sistemas altamente no lineales. La adición del Jacobiano también puede hacer que el modelo sea más difícil de diseñar y, como tal, a pesar de sus ventajas, puede ser más apropiado implementar el KF por simplicidad y interoperabilidad.

A menos que se indique lo contrario, todas las imágenes son del autor.

We will continue to update Zepes; if you have any questions or suggestions, please contact us!

Share:

Was this article helpful?

93 out of 132 found this helpful

Discover more

Inteligencia Artificial

Investigadores de UCSD y Microsoft presentan ColDeco una herramienta de inspección sin código para columnas calculadas.

En el artículo “COLDECO: una herramienta de inspección de hojas de cálculo para código generado por IA” u...

Inteligencia Artificial

Cómo gané en el fútbol de fantasía italiano ⚽ utilizando el aprendizaje automático

Como ingeniero mecánico con un gran interés en la programación y la informática, me fasciné por el mundo del aprendiz...