Recordando al Filtro Kalman

Resumen

En 1960 R.E. Kalman publicó su famoso documento donde describía una solución recursiva al problema de filtrado de datos discretos. Desde ese entonces, el filtro de Kalman fue sujeto de muchísimas investigaciones y aplicaciones particularmente en el área de la navegación autónoma o asistida. Para una mayoría de los mortales es prácticamente imposible comprender el significado completo de lo que representa y significa el filtro Kalman a partir de definiciones y complicadas ecuaciones. En el presente trabajo se muestra las ecuaciones que representan a este filtro además de su aplicación al caso de un motor de corriente directa empleando para la simulación el software MATLAB.

Palabras clave: Filtro Kalman, Estimación de estados, Estimación óptima.

Introducción y principios de la estimación de estados

El 19 de Mayo, es el aniversario del nacimiento de éste Ingeniero Húngaro-Americano. Matemático teórico y ex profesor Universitario congratulado por el presidente Norteamericano Barack Obama el 7 de Octubre del 2009 por sus notables aportes y la co-invención del filtro Kalman, una formulación matemática muy empleada en sistemas de control, aviónica, astronomía e incluso en economía. Las ideas de Kalman sobre el filtro fueron inicialmente recibidas con escepticismo, tanto que fue forzado a publicar sus resultados en una revista para ingenieros mecánicos (ASME). Va como homenaje a este Ingeniero, el presente trabajo.

El primer método para realizar estimaciones óptimas a un conjunto de datos con ruido /1/, es el método de los mínimos cuadrados. Descubierto por Friedrich Gauss en 1795. Este método fue muy empleado para problemas de estimación lineal en astronomía.

En 1801 El astrónomo Italiano Giuseppe Piazzi estuvo observando un catálogo de estrellas. Sin saberlo, Piazzi incluyo un error de impresión. Mientras buscaba una estrella perdida, Piazzi descubrió en su lugar algo que lo movió posteriormente. Fue el planeta enano Ceres, el cuerpo más grande en la faja de asteroides y el primero a ser descubierto, pero Piazzi no conocía aquello aún. El fue capaz de seguirle la pista y medir su movimiento aparente en referencia a una estrella "fija" de fondo durante 41 noches antes de que se haya acercado al sol y desaparecido.

El 24 de enero de 1801, Piazzi escribió sus descubrimientos y los envió a Johann Bode, mejor conocido por su denominada ley de Bode, la cual establece que las distancias de los planetas desde el sol en unidades astronómicas (una unidad astronómica es igual a la distancia entre la Tierra y la Luna) está dada por la secuencia:

Actualmente se conoce que no fue Bode, sino Johann Tietz quien primero propuso esta fórmula en 1772. En aquel entonces solo se conocía de la existencia de seis planetas. En 1781, Friederich Herschel descubrió Urano, el cual encaja cercanamente en esta fórmula para n=6. No se descubrió planeta alguno para n=3. Estimulados por Bode, una asociación de astrónomos Europeos estuvo buscando el octavo planeta "perdido" por cerca de 30 años. Piazzi no fue parte de esta asociación.

La carta de Piazzi a Bode, no llegó sino hasta el 20 de marzo (los e-mail fueron descubiertos muchísimos años después). Bode sospechó que el descubrimiento de Piazzi podría estar relacionado con el planeta perdido pero existían datos insuficientes para determinar sus elementos orbitales por los métodos entonces disponibles. Es un problema donde participan ecuaciones no lineales que el mismo Newton declaró ser uno de los más complicados en astronomía matemática. Nadie resolvió el problema, como resultado, Ceres se perdió en el espacio nuevamente.

Los descubrimientos de Piazzi no fueron publicados hasta el otoño de 1801. El posible descubrimiento – y subsecuentes pérdidas – de un nuevo planeta coincidió con el inicio de un nuevo siglo, fueron novedades excitantes. Contradecían a una justificación filosófica de solo siete planetas cuyo número se conocía antes de Ceres y defendida por el respetable filósofo Georg Hegel, en medio de otros. Hegel había publicado recientemente un libro en el cual el castigaba a los astrónomos por enfermar esos tiempos buscando un octavo planeta cuando en ese entonces dominaban las justificaciones filosóficas de que sean solo siete. El nuevo objeto celestial fue sujeto de conversaciones en círculos intelectuales por todo lado. Afortunadamente, el problema llamo la atención de un matemático Gótico de 24 años, llamado Carl Friedrich Gauss.

Gauss dedicó mucho de su tiempo al problema de la determinación de la órbita de Ceres en Diciembre y envió sus resultados a Piazzi. El nuevo planeta (luego reclasificado como asteroide el 2006) fue encontrado en el último día del año.

Gauss no publicó sus métodos de determinación de la órbita hasta 1809. En esa publicación, el también describe el método de los mínimos cuadrados que el había descubierto en 1795 a la edad de 18 años y que fue utilizada refinadamente para la determinar la órbita de Ceres.

El método de los mínimos cuadrados fue sujeto de continuo interés y beneficios para generaciones de científicos y tecnólogos. Fue el primer método de estimación óptima y probó una importante conexión entre las ciencias teóricas y experimentales: Proporcionó a los experimentalistas un método práctico para estimar los parámetros desconocidos de modelos teóricos.

Algunas definiciones

Filtro: En términos generales un filtro es aquello que discrimina elementos indeseables de un conjunto. En sistemas de control, un filtro es aquel que discrimina determinadas frecuencias de una señal eléctrica que pasa por él. Ejemplo oleo-hidráulico: el aceite de un motor que atraviesa por el filtro de aceite el cual detiene las partículas en suspensión con un tamaño superior a los microscópicos orificios de dicho filtro.

Proceso: conjunto de actividades que se realizan bajo ciertas circunstancias y con un fin determinado.

Ruído: Es una señal sin información ni datos y que perturba a otras señales que llevan información o datos. Ejemplo la señal del sonido en el auricular de un teléfono perturbado por el ruido de un motor de electrodoméstico.

Ruído Blanco: Una señal aleatoria cuya densidad espectral de potencia es una constante. Los valores temporales en dos tiempos distintos no guardan correlación estadística.

Señal: Dentro el campo que nos compete, una señal es la variación en el tiempo de alguna magnitud física y que se emplea para transmitir información.

Sistema: para el presente caso, es cualquier dispositivo destinado a la generación, transmisión, procesamiento o almacenamiento de señales.

Definición del problema

Suponga que tiene una señal, por ejemplo la señal de voz proveniente de un micrófono, la señal proveniente de un radar que contiene la posición de un objeto en el aire, la señal proveniente de los electrodos de un equipo de ECG, etcétera. Y por supuesto estas señales están contaminadas con ruido proveniente de su entorno. Como extraer este ruido?. La respuesta es hacer uso de filtros. Antiguamente se empleaban filtros analógicos, realizados con dispositivos tales como capacitores, resistencias, inductores y amplificadores operacionales básicamente; sin embargo hoy en día y gracias a las modernas técnicas de adquisición y procesamiento de señales, este proceso se realiza totalmente con sistemas microprocesados. Luego una vez adquirida la señal en forma digital, queda el procesarla digitalmente para reducir a prácticamente cero el ruido contenido en la señal que previamente adquirimos o medimos.

Una vez obtenida la señal discreta con ruido, la idea rápida que podría venir a mente es tomar el promedio de las muestras consecutivas, pero esta aproximación no trabaja en la mayoría de los problemas de la vida real.

Existe varias técnicas para el Procesamiento Digital de Señal (PDS), una de las cuales es la que nos ocupa en el presente: el filtro Kalman.

Sistemas lineales

La mayoría de los procesos conocidos por el hombre, pueden ser descritos mediante un conjunto de ecuaciones lineales. Para el caso que nos ocupa, un sistema lineal es simplemente un proceso que puede ser descrito por las siguientes ecuaciones:

A, B, H son matrices, k es el subíndice que señala un evento en el tiempo discreto, x es llamado vector de estados del sistema, u es un vector de entradas conocidas al sistema, z es el vector de salidas o medidas en el sistema, w y v son vectores de ruido asociados al proceso y a las medidas respectivamente. La ecuación (1) es mas conocida como ecuación de estado y la (2) es la ecuación de la salida. El vector x contiene toda la información relacionada con el estado del sistema pero no puede ser medido directamente. En su lugar se mide z el cual es una función de x que está contaminado con ruido z.

La teoría y algoritmo del filtro Kalman

Conocido el sistema lineal que describe el proceso, nuestro objetivo ahora es estimar el estado x del sistema a partir de las mediciones disponibles z. Como podemos obtener la mejor estimación del estado x?. Lo deseable es obtener una estimación precisa del verdadero estado teniendo en mente que no podemos medir directamente los mismos.

Que criterios deberían satisfacerse?. Primero nuestro filtro debe hacer que el valor promedio del estado estimado sea igual al valor promedio del verdadero estado; esto matemáticamente significa que el valor esperado de la estimación debe ser igual al valor esperado del estado. Segundo, nuestro filtro debe hacer que el estado estimado varíe poquísimo del estado verdadero; esto significa matemáticamente que deseamos un estimador con la más pequeña varianza en el error.

Sucede que el filtro Kalman es el estimador que satisface esos dos criterios; sin embargo el filtro Kalman no puede aplicarse a menos que se satisfagan ciertas condiciones en el ruido que afecta a la señal. La primera condición es que el valor promedio de las señales de ruido w y v sea cero y la segunda es que no hay correlación entre dichas señales; es decir para algún instante k, las variables aleatorias w y v son independientes. Se definen las matrices de covarianza de ruido Q y R como:

Donde T significa la traspuesta de un vector o matriz y E(*) simboliza el valor expectado o esperado (esperanza matemática).

Existen varias formas de expresar las ecuaciones del filtro Kalman, una de ellas es /3/:

En estas ecuaciones, el superíndice (-1) significa inversión matricial, un (T) significa traspuesta. La matriz K se denomina ganancia Kalman y la matriz P es la covarianza del error de estimación, el estado estimado es (algunos suelen denominar "sombrero" al circunflejo).

La ecuación (8) es intuitiva y consta de la suma de dos términos. El primer término de la suma es empleado para hallar el estado estimado en el instante k; el segundo término de la suma no es más que un término corrector que nos indica en cuanto se debe corregir dicho estado estimado debido a nuestras mediciones.

Las ecuaciones 5 y 6 son también conocidas como ecuaciones de actualización en tiempo o ecuaciones de predicción; las ecuaciones 7, 8 y 9 son conocidas como ecuaciones de actualización de mediciones o ecuaciones correctoras. Para iniciar el proceso, necesitamos conocer inicialmente x0 y P0. Luego se puede ejecutar el algoritmo de acuerdo a lo que señala la figura 1; los datos iniciales ingresan a un "circuito cerrado" hasta que uno decida detenerlo.

Una vez detenido el proceso de cálculo, se obtiene la estimación almacenada de las anteriores iteraciones.

La parte mas penosa es determinar las matrices R y Q. Generalmente R es fácil de encontrar debido a que podemos medir el proceso de alguna forma (mientras está trabajando el filtro) de tal forma que podamos tomar algunas mediciones "off-line" para así determinar la varianza de ruido en las mediciones. Pero encontrar Q no es tan obvio pues típicamente no tenemos la posibilidad de observar directamente el proceso que estamos estimando.

Ejemplo de aplicación

Para probar el filtro, rescatamos el modelo de segundo orden de un motor de corriente directa presentado en /6/. El modelo continuo es:

La figura 2, muestra la evolución de la velocidad del rotor en el tiempo, para un período de análisis de 0.1 s.

Para emplear el filtro, asumiéremos lo siguiente: se trata de un motor DC donde la corriente fue medida con un error de 0.2 A (desviación estándar), la tensión de entrada es de 12 V con un ruido de 1.1 V (desviación estándar).

De modo a obtener las matrices Q y R, podemos plantear la siguiente aproximación. Como la desviación estándar del ruido en las mediciones es de 0.2, la matriz R será simplemente R = (0.2)2 = 0.04. Ahora necesitamos obtener Q de forma aproximada, para ello diremos que como la corriente de armadura es proporcional a 0.2802 veces la tensión aplicada al motor y el ruido en la tensión aplicada es 1.1 V, entonces la varianza del ruido en la corriente de armadura será aproximadamente (0.2802)2(1.1)2 = 0.095. Similarmente, puesto que la velocidad es proporcional a 0.3521 veces la tensión aplicada, la varianza en el ruido de la velocidad será (0.3521)2(1.1)2 = 0.15. Finalmente la covarianza del ruido en la corriente de armadura y ruido en la velocidad será igual a la desviación estándar del ruido en la corriente de armadura multiplicada por la desviación estándar del ruido en la velocidad del rotor, con esto (0.2802*1.1)(0.3521*1.1) = 0.1194. Luego en forma matricial, podremos escribir:

El anexo, posee un programa en MATLAB para realizar la simulación correspondiente. Observe que hay un "factor", el cual inicialmente es 1, pero que con ese valor la simulación no es satisfactoria. Si se cambia dicho factor a 5, 10, 100, podemos ver en las figuras 3, 4, 5 y 6 como mejora el estado estimado, el cual se aproxima a la solución presentada en la figura 2.

El lector puede observar que cuanto mas grande es el "factor", la matriz Q es más "pequeña" y consecuentemente la estimación del estado es más satisfactoria aproximándose al resultado "ideal" de la figura 2. Se deja al lector intentar el caso para cuando en lugar de hacer mas "pequeño" a Q, se hace más "grande" a R.

Conclusiones

A partir de mediciones con ruido y teniendo el modelo del proceso, hemos estimado el estado del mismo, que en este caso es la velocidad del rotor.

Cabe acá hacer las siguientes puntualizaciones:

• El filtro Kalman no es del todo un filtro, sino un estimador.

• Es un método recursivo, fácilmente implementable como algoritmo computacional pues emplea sus resultados previos como datos de entrada para nuevos resultados.

• El filtro Kalman asume un sistema dinámico lineal y uno observable corrompido con ruido blanco.

• Permite la estimación de estados presentes, pasados y futuros en condiciones de incertidumbre.

Anexo

Programa en MATLAB para realizar la simulación (también se puede copiar el script en http://electricosas.blogspot.com/).

%DATOS DEL SISTEMA

A=[0.4146 -0.0066;1.4643 0.9916];

B=[0.2802;0.3521];

H=[1 0];

u=12;

%MATRICES DE VARIANZA Y ERRORES

%cambiar este factor por 5, 10, 100 y ver

%como mejora la velocidad estimada

%—————-

factor=1;

%—————-

Q=[0.095 0.1194;0.1194 0.15]/factor;

R=0.04;

err_med_corr_armad=0.2;

err_med_tens_ent=1.1;

%DATOS INICIALES

x_est_apriori=[0;0];

P_apriori=[0 0;0 0];

K=[0;0];

x_est=[0;0];

P=Q;

x=[0;0];

%PARA LA SIMULACION

ia=zeros(1,1000);

w=zeros(1,1000);

t=1:1000;

t=t*0.0001;

%PROCESO ITERATIVO

for i=1:1000

%GENERAMOS LAS MEDICIONES RUIDOSAS

%—————————————————————–

ruido_de_proceso=err_med_tens_ent*[0.2802*randn;0.3521*randn];

x=A*x+B*u+ruido_de_proceso;

ruido_mediciones=err_med_corr_armad*randn;

z=H*x+ruido_mediciones;

%—————————————————————–

%ECUACIONES DE PREDICCION

%—————————————————————–

x_est_apriori=A*x_est+B*u;

P_apriori=A*P*A'+Q;

%—————————————————————–

%ECUACIONES DE CORRECCION

%—————————————————————–

INNOVA=H*P_apriori*H'+R;

K=P_apriori*H'/INNOVA;

x_est=x_est_apriori+K*(z-H*x_est_apriori);

P=A*P_apriori*A'-A*P_apriori*H'*(1/R)*H*P_apriori*A'+Q;

%—————————————————————–

%ALMACENAMOS RESULTADOS DE CADA ITERACION

%—————————————————————–

ia(i)=x_est(1);

w(i)=x_est(2);

%—————————————————————–

end

%GRAFICAMOS LA VELOCIDAD DEL ROTOR

plot(t,w,'k')

Bibliografía

/1/ Kalman Filtering, Theory and Practice Using MATLAB, Third Edition, Mohinder S. Grewal, Angus P. Andrews, John Wiley and Sons Inc., 2008, ISBN 978-0-470-17366-4.

/2/ Kalman Filter, Vedran Kordic, Intech, May 2010, ISBN: 9533070940.

/3/ Optimal State Estimation, Dan Simon, John Wiley and Sons Inc., 2006, ISBN -100-471-7085 8-5.

/4/ A Kalman Filter Primer, R. L. Eubank, Chapman and Hall/CRC, 2006, ISBN -100-8247-2365-1

/5/ Kalman filter and Neural Networks, Simon Haykin (editor), John Wiley and sons, 2001, ISBN 0-471-22154-6.

/6/ A Tutorial on Dynamic Simulation of DC Motor and Implementation of Kalman Filter on a Floating Point DSP, Padmakumar S., Vivek Agarwal, and Kallol Roy, Mumbai, India, World Academy of Science, Engineering and Technology 53, 2009.

BIOGRAFIA.

Rogelio José Choque Castro, Nacido en La Paz Bolivia, es Ingeniero Electricista titulado en la Universidad Mayor de San Andrés (UMSA). Trabajó durante diez años en la Industria Textil (área de mantenimiento), Supervisor Proyecto de Electrificación Rural Illimani Sud, Residente de Obra Proyecto Porvenir-Chive, Residente y Superintendente para la empresa ENERGO S.R.L., Supervisor proyecto MCH-Totorapampa, Supervisor Queñi-Cotosi, fue docente de la Universidad de Aquino Bolivia (UDABOL) y Universidad Los Andes. Sus áreas de interés: Simulación de Transitorios en Sistemas de Potencia, Electrónica de Potencia, Procesamiento Digital de Señales, Sistemas de Control.

Página: http://electricosas.blogspot.com/

Autor:

Rogelio José Choque Castro