El Filtro de KalmanAntonio Miguel Fumero Reverón

Temas Especiales de Análisis NuméricoJosé Manuel Corrales Sendino

Madrid. Día 18, Febrero de 2000

________________________________________________________________TEAN

Motivación y Objetivo

La elaboración de este documento parte del interés personal del alumno por ciertos

temas relacionados con las comunicaciones, y los problemas clásicos que ha planteado

la comunicación entre entidades (máquinas) dotadas de cierta inteligencia (entendida

en el sentido de la posesión de una capacidad de decisión determinada). Este interés

surge en cierta medida tras la lectura de los artículos de Shannon acerca de la teoría

matemática de la comunicación, y sobre todo, de la conocida obra de Norbert Wiener

“Cibernética”, originalmente editada en 1948 por el MIT bajo el título de “Cybernetics or

Control and Communication in the Animal and the Machine”. Llama la atención la

analogía entre los problemas de la predicción estadística y el control determinístico; así

como la cantidad de trabajo que ha surgido en torno a las ecuaciones normales de

Wiener.

El interés por el tema se ha visto incrementado tras la asistencia a un pequeño

seminario acerca de los sistemas no lineales y de control adaptativo que se impartió por

vez primera como asignatura de libre elección en el P94 durante el segundo

cuatrimestre del curso 1998/1999. En el seminario se presentaba el problema de

Wiener y sus aplicaciones en el mundo del control adaptativo y ecualización de canales

de comunicación.

El objetivo principal del documento es volver a plantearse esos problemas, así como las

soluciones que propuso en su día R. E. Kalman, que han dado lugar a una serie de

algoritmos numéricos de amplia aplicación en diversos campos de nuestra técnica. Para

ello se ha acudido a algún texto de filtrado adaptativo y ecualización, así como a los

artículos originales de Kalman donde propone su solución del problema de Wiener y

sus posibles aplicaciones.

Se han pretendido usar las ideas y conceptos presentados en el curso de TEAN,

impartido en el primer semestre del presente curso 1999/2000, para elaborar un estudio

del problema y las aplicaciones del Filtro de Kalman con el objetivo de que sirva de

______________________________________________________________________

________________________________________________________________TEAN

alguna forma para la evaluación de la asignatura en términos de su mayor o menor

grado de aprovechamiento.

Dicho sea ya antes de nada que el autor se declara una persona más preocupada por

los conceptos involucrados en el problema que se presenta y su alcance en cuanto a

sus implicaciones en las técnicas de uso común en la ingeniería, que por el detalle

operativo de los métodos numéricos que se emplean en la resolución de los problemas

concretos. Este hecho, quizás motivado por la lejanía de los primeros años de

universidad en que uno se acostumbraba a lidiar en el día a día con el utillaje y la

operativa matemática, podría llevar a pensar que estas líneas no constituyen otra cosa

que un de ejercicio de demagogia o disquisición metafísica sin utilidad práctica alguna.

______________________________________________________________________

________________________________________________________________TEAN

Introducción

Aunque la difusión que va a recibir este documento no va a ser en absoluto amplia, creo

que no es mala idea comenzar con una pequeña introducción al problema del filtrado

adaptativo, estimación lineal, y optimización que servirán al menos para establecer

unos términos básicos sobre los que se desarrollará todo el documento.

Lo primero es dejar claro que cuando hablamos de filtros adaptativos estamos hablando

de sistemas lineales programables; es decir de filtros lineales que operan sobre una

serie temporal de entrada, pero cuyos coeficientes (que definen la operación de filtrado

que realizan sobre esa entrada) se pueden variar a voluntad de un instante temporal al

siguiente (de una muestra a la siguiente) mediante un algoritmo adaptativo que se

basará en alguna relación recursiva:

hi(n+k) = hi(n) + g{x(n),ŷ(n),e(n)}

donde hi(n) es el coeficiente i-ésimo del filtro en el instante temporal n, x(n) es un vector

de las entradas pasadas del filtro, ŷ(n) es un vector de las salidas pasadas del filtro, y

e(n) = y(n) - ŷ(n) es un vector que representa el pasado de la señal de error.

La estructura básica que se representa en la Figura 1 tiene tres modos de

funcionamiento fundamentales:

______________________________________________________________________

Σ- +

Filtro Programable

Algoritmo Adaptativo

h(n)

ŷ(n)

y(n)

x(n)

e(n)

Figura 1:Diagrama de Bloques de un Filtro Adaptativo Genérico (Tomada de [Mulgrew88])

________________________________________________________________TEAN

(i) Identificación de Sistemas (Figura 2(i)) : se trata de alimentar al sistema

desconocido que se pretende identificar, y a la estructura adaptativa de

antes con la misma serie temporal de entrada; de forma que tras la

convergencia la salida del filtro adaptativo ŷ(n) se aproximará a la salida

del sistema desconocido y(n) en un sentido óptimo (normalmente se usará

el criterio de convergencia en media cuadrática). En esa situación se dirá

que hemos identificado el sistema. Normalmente tendremos ruido blanco

en la entrada (errores en las observaciones) que nos alejará del

comportamiento ideal que se pretende.

(ii) Problema Inverso de Análisis o Modelado de Sistemas(Figura 2(ii)) : es

típico de las estructuras que se usan para la ecualización de canales de

comunicación. En este caso se trata de que la salida del sistema

desconocido sirva de entrada para el filtro adaptativo. Lo que se pretende

es obtener la salida del filtro adaptativo una versión retardada de la señal

de entrada, con lo cual, tras la convergencia, la función de transferencia

del filtro coincidirá con la inversa de la función de transferencia del sistema

desconocido que lo precede (típicamente el canal de comunicaciones).

(iii) Predicción Lineal(Figura 2(iii)) : se usa mucho en los codificadores en

canales de comunicación vocal. En este caso le metemos como entrada al

filtro adaptativo una versión retardada de la misma señal que deseamos

tener a la salida del mismo; con lo cual se trata de predecir las muestras

futuras de la entrada. Este comportamiento sólo se consigue si la señal de

entrada está bastante lejos de tener un espectro plano.

Con estas ideas más o menos claras y algunos conceptos de sistemas lineales que se

han ido destilando a lo largo de los años pasados en la Escuela, pienso que sería

buena idea plantear el problema de la estimación lineal óptima y la ecuación de Wiener

que es la motivación principal de toda esta historia.

El problema de la estimación lineal se puede expresar de una forma más o menos

rigurosa con un enunciado como este: Dada una secuencia aleatoria de observaciones

______________________________________________________________________

________________________________________________________________TEAN

{x(n)}, que se corresponde con una versión distorsionada de una señal {y(n)} (otra

secuencia aleatoria que transporta cierta información), encontrar un filtro lineal que

opere sobre {x(n)} para obtener una estimación { ŷ(n)} de {y(n)}. La calidad de la

estimación se medirá con una función f(.) de la secuencia de error {e(n)}, considerado

este como la diferencia e(n) = y(n) - ŷ(n), que asignará un precio o penalización a las

estimaciones incorrectas (de ahí que habitualmente se la conozca como función de

pérdidas). Se trata de una función que debe ser positiva (o no negativa), y decreciente.

Puesto que la secuencia de error también es una secuencia aleatoria, se suele elegir el

filtro de forma que minimice una función de coste l(.) que sea la esperanza matemática

de la función de pérdidas: l(e(n)) = E[f(e(n))] . La función de coste más habitual suele

ser el error cuadrático medio (MSE) (n) = E[e2(n)].

El filtro FIR óptimo se obtendría siguiendo estas consideraciones:

ŷ(n) es la salida de un filtro causal, y como tal se puede expresar como la suma de

convolución de la secuencia de entrada y la respuesta al impulso, que se puede truncar

para obtener un filtro FIR de orden N-1

ŷ(n) = hix(n-i) ; (0iN-1)

que si se expresa en forma vectorial es algo como

ŷ(n) = hT x(n) ; h = [h0 h1 ... hN-1]T ∧ x(n) = [x(n) x(n-1) ... x(n-N+1)]T

La función de coste MSE la podemos poner como

= E[y2(n)] + hTxxh – 2hTxy

xx = E[x(n)xT(n)] ℳNN es la matriz de autocorrelación

______________________________________________________________________

Sistema Desconocido

Filtro AdaptativoΣ

- +

x(n) y(n)

Figura 2(i):Identificación de Sistemas (Tomada de [Mulgrew88])

________________________________________________________________TEAN

xy = E[x(n)y(n)] ℳN1 es el vector de correlación cruzada

Para minimizar la función de coste, igualaremos el gradiente a cero:

= ∂/∂h = 2xxh – 2xy = 0 ⇨ xxhopt = xy

Resulta que si la densidad espectral de potencia de la secuencia de entrada {x(n)} no

tiene nulos, la matriz de autocorrelación es definida positiva, y por tanto no singular. Y

sabemos que bajo esas condiciones la respuesta al impulso óptima es única y viene

dada por

hopt = -1xx xy⇨opt = E[y2(n)] - hT

optxy

Esa ecuación es la que se conoce como el filtro FIR de Wiener (o filtro de Levinson). Lo

que ocurre aquí es que necesitamos la matriz de autocorrelación y el vector de

correlación cruzada; es decir necesitamos estadísticos de segundo orden cuando lo que

normalmente se conoce son las secuencias aleatorias de datos. La labor de determinar

el filtro óptimo a partir de esos datos se le suele encomendar a un filtro adaptativo, que

se puede ver en esos términos como un algoritmo que opera sobre las secuencias

{x(n)} e {y(n)} para formar a partir de ellas un vector de respuesta al impulso h(k)

variable con el tiempo que converge en media a la respuesta al impulso óptima. El filtro

de Levinson será por tanto el objetivo que debe perseguir el filtrado adaptativo.

El cálculo de la hopt requiere la solución de N ecuaciones lineales simultáneas con N

incógnitas. Para una matriz no singular en general el método más eficiente sabemos

que es el de eliminación de Gauss, que requiere O(N3) operaciones, siendo también N

el número de coeficientes del vector de respuesta al impulso. De lo que se dio cuenta

Levinson fue del hecho de que la matriz de autocorrelación tiene una estructura muy

especial:

(i) es simétrica ⇨ xx[i,j] = xx[j,i]

______________________________________________________________________

Sistema Desconocido

Filtro AdaptativoΣ

- +x(n)

Figura 2(ii):Problema Inverso de Modelado o Análisis de Sistemas (Tomada de [Mulgrew88])

________________________________________________________________TEAN

(ii) y es de Toeplitz⇨ xx[i,j] = xx[i-j]

siendo xx[i,j] el elemento en la fila i, columna j de la matriz de autocorrelación, y {xx(m)}

la secuencia de autocorrelación asociada con la secuencia estacionaria en sentido

amplio {x(n)}. Levinson se dio cuenta de que con esta estructura sólo necesitaba la

primera fila para definir toda la matriz; con lo cual este hombre ideó un algoritmo que

sólo requería O(N2) operaciones.

Hay también otra forma más intuitiva de atacar el problema que se conoce como SMI

(de Sampled Matrix Inversion) que consta de dos fases: en una primera fase se hace

una estimación de los coeficientes de correlación y autocorrelación a partir de las

secuencias de datos disponibles; y en la segunda fase se usan esos coeficientes

estimados para construir la ecuación de Wiener, que se resuelve luego con el algoritmo

de Levinson.

Cuando sólo tenemos conjuntos finitos de datos (x(n), y(n), n=0,...k), se sustituye el

operador esperanza matemática de la función de coste MSE por un sumatorio,

obteniéndose lo que se llama la función de coste LS ( por Least Squares). Al final

resulta que el problema de la minimización LS es dual del de la minimización MSE

quedando algo como

rxx(k)h(k) = rxy(k)

rxx (k) = (n=0,...,k)x(n)xT(n) ℳNN es la matriz de autocorrelación

rxy (k) = (n=0,...,k)x(n)y(n) ℳN1 es el vector de correlación cruzada

Existirá solución única si existe rxx (k) y es no singular.

Para kN-1, rxx (k) será siempre singular.

Para kN, la singularidad de rxx (k) dependerá de los conjuntos de datos particulares

que tengamos disponibles.

El problema de la singularidad aveces se puede soslayar usando la pseudoinversa de

Moore-Penrose, y hay autores que lo han propuesto.

Hay varias formas de construir estimaciones LS dependiendo de las hipótesis que se

hagan sobre los datos disponibles: se puede usar directamente la covarianza, la

autocorrelación, o el enventanado (pre- o post- enventanado).

También existe un algoritmo recursivo RLS para los casos en que se hace necesario

actualizar la estimación LS según vamos teniendo nuevos datos disponibles.

______________________________________________________________________

________________________________________________________________TEAN

Hay toda una serie de algoritmos, que se conocen como algoritmos rápidos. Surgieron

a partir de la publicación del FKA (Fast Kalman Algorithm) que trataba de aprovechar la

estructura peculiar de la matriz rxx (k) de la estimación LS; le siguieron más algoritmos

buscando una mejora en la eficiencia computacional: FAEST (Fast A poSteriori Error

Technique), FTF (Fast Transversal Filter).

También los métodos estocásticos de gradiente dan técnicas iterativas donde el número

de operaciones que se deben llevar a cabo para encontrar la solución no se conoce de

antemano. El método de descenso en el sentido del gradiente es una técnica iterativa

que se usa en programación lineal y problemas de optimización para encontrar una

variable que minimice una función de coste objetivo. Para aplicarlo a la ecuación de

Wiener: se hace una estimación de la solución hi (siendo i el paso en el proceso de

iteración), que tendrá asociada un valor de la función de coste MSE, (hi). Luego se

hace una mejora de la estimación en dos fases: en la primera se calcula el gradiente de

la misma, i que será un vector que irá en la dirección de la máxima pendiente de la

función de coste MSE; en la segunda fase se construye una nueva estimación

restándole a la primera una versión escalada del gradiente (hi). Queda algo así:

hi+1 = hi - (hi)

i = ∂/∂h(hi) = 2xxhi – 2xy

como condición de convergencia tenemos que

0<<1/máx

siendo máx el mayor autovalor de la matriz de autocorrelación

Este método evita la necesidad de la inversión de la matriz de autocorrelación, pero

sigue necesitando conocer esa matriz y el vector de correlación cruzada. Pero hay un

algoritmo que no necesita el conocimiento de esos estadísticos, es el LMS (por Least

______________________________________________________________________

Filtro AdaptativoRetardo Σ

- +

x(n)

Figura 2(iii):Predicción Lineal (Tomada de [Mulgrew88])

________________________________________________________________TEAN

Mean Squares). Este algoritmo sustituye el proceso de iteración por una recursión

temporal y el gradiente por una estimación del mismo.

h(k+1) = h(k) - (k)

(k) = ∂/∂h(h(k)) = ∂/∂h E[(y(n) – hT(k)x(n))2] = -2E[x(n)(y(n) – hT(k)x(n))]

(k) = -2x(k+1)e(k+1)

e(k+1) = y(k+1) – hT(k)x(k+1)

Al final el algoritmo se puede expresar como

h(k+1) = h(k) + 2x(k+1)e(k+1)

que es bastante más simple que el RLS.

Existe una variante, el BLMS, en donde se mantiene constante la estimación de la

respuesta al impulso en un bloque de L puntos de la secuencia de datos. Da muy

buenos resultados cuando se usa junto con la FFT.

El Filtro de Kalman

Una vez que estamos de acuerdo sobre donde comienza toda esta historia, si nos

vamos ahora a mirar en cualquier texto dedicado al filtrado adaptativo y sus

aplicaciones, y buscamos alguna referencia al filtro de Kalman, encontraremos en

______________________________________________________________________

________________________________________________________________TEAN

muchos su utilización como Ecualizador IIR adaptativo justificada más o menos en los

siguientes términos.

El filtro IIR de Wiener, que resulta ser bastante mejor que el FIR para canales

conocidos, comienza a dar problemas cuando, ante un canal desconocido o variable en

el tiempo exige la utilización de estructuras adaptativas. Si la estimación se realiza con

el mismo algoritmo LMS que se usa para el caso FIR, existe cierta probabilidad de que

los polos del filtro se salgan de la circunferencia unidad; pudiendo dar lugar a una

situación de inestabilidad si se mantienen así durante mucho tiempo.

La idea es que la realización del ecualizador IIR óptimo de Wiener implica

consideraciones y conceptos relacionados con el espacio de estados, y eso es lo que

nos lleva al planteamiento original que Kalman presenta en sus artículos de principios

de los años 60.

Lo que se hace es usar el Filtro de Kalman como ecualizador; y para darle carácter

adaptativo a la estructura se usa en paralelo con un algoritmo de identificación de

sistemas que nos dará los M coeficientes de la respuesta al impulso del canal (una

estimación que se minimizará con un criterio MSE). También hay autores que han

propuesto soluciones donde meten los coeficientes del canal en la información asociada

al estado del modelo de Kalman (es lo que llaman el Filtro de Kalman Extendido, EKF).

El modelo de Kalman se puede poner, en tiempo discreto y con una notación más o

menos estándar como sigue.

S(k+1) = A(k+1)s(k) + w(k+1)

z(k) = H(k)s(k) + u(k)

s(k) es lo que se llama el vector de estado y contiene los valores de los M parámetros

que definen el estado del sistema en el instante k.

A(k) es una matriz de dimensiones MM, que se conoce como la matriz de transición

de estado.

______________________________________________________________________

________________________________________________________________TEAN

H(k) es una matriz de dimensiones LM, que se conoce como la matriz de

observaciones.

w y u son dos vectores de M y L componentes respectivamente; y representan a dos

procesos blancos incorrelados y de media nula de covarianzas W y U, de forma que:

E[w(k)wT(k)] = W(k)

E[u(k)uT(k)] = U(k)

Se supone que las matrices A, W, H, U, aunque pueden variar con el tiempo, son

conocidas a priori.

Las estimaciones óptimas del vector de estado s(k) se generan de forma recursiva a

partir de la secuencia de observaciones ruidosas {z(k)}. Las ecuaciones pueden

ponerse como sigue:

ESTIMACIÓN

Ŝ(k/k) = Ŝ(k/k-1) + K(k)[z(k) – H(k) Ŝ(k/k-1)]

PREDICCIÓN

Ŝ(k/k-1) = A(k) Ŝ(k-1/k-1)

GANANCIA

K(k) = V(k/k-1)HT(k)[H(k)V(k/k-1)HT(k) + U(k)]-1

COVARIANZA DEL ERROR

V(k/k-1) = A(k)V(k-1/k-1)AT(k) + W(k)

V(k-1/k-1) = V(k-1/k-2) – K(k-1)H(k-1)V(k-1/k-2)

V(k/k) = E[(s(k) – Ŝ(k/k))(s(k) – Ŝ(k/k))T]

V(k/k-1) = E[(s(k) – Ŝ(k/k-1))(s(k) – Ŝ(k/k-1))T]

Al final, ¿Qué es lo que pasa aquí? ¿Qué es lo más molesto en el filtro IIR de Wiener?

Conseguir la factorización de fase mínima del espectro de potencia de las

observaciones.

______________________________________________________________________

________________________________________________________________TEAN

Si todos los procesos son estacionarios, y el ruido en las observaciones es blanco, el

filtro de Kalman para el estado estacionario y el filtro IIR de Wiener son idénticos.

Un filtro FIR nos lleva directamente a una representación sobre la base del espacio de

estados. Pero hay que tener cuidado: un vector de estado con (M-1) estados nos llevará

a una formulación donde el ruido de planta y el de las observaciones están correlados,

y un filtro de Kalman en esa situación es condicionalmente estable. Con M estados

representando al canal, tendríamos sin embargo una formulación en la que el ruido de

planta y el de las observaciones estarían incorrelados, y el filtro de Kalman sería

incondicionalmente estable.

Para montar un ecualizador IIR con el filtro de Kalman tendríamos la siguiente

formulación:

s(k) = [s(k) s(k-1) ... s(k-M+1)]T es el estado del modelo FIR de canal

s(k) = as(k-1) + bs(k)

a es una matriz de desplazamiento de MM

bT = [1 0 0 ... 0]

↪x(k) = hT(k)s(k) + n(k)

{s(k)} y {n(k)} son secuencias de ruido blanco, incorreladas, y de media nula con

varianzas s2(k) y n

2(k).

Este es un planteamiento válido para cuando sólo necesitamos una estimación del

estado actual del canal. Pero si tratamos con canales que no son de fase mínima,

necesitaremos montarnos otra historia:

s(k) = [s T (k) s T (k-1) ... s T (k-d)]T

El vector de estado se segmenta para contener (d+1) elementos

s T (k) = [s(k) s(k-1) ... s(k-M+1) ... s(k-d)]T

Los vectores a y b se construyen igual que antes, pero con (d+1) elementos

La matriz de observaciones se construye rellenando con ceros la respuesta al impulso

del canal:

H(k) = [h0(k) h1(k) ... hM-1(k) 0 0 ... 0]

______________________________________________________________________

________________________________________________________________TEAN

Con estas consideraciones se pueden escribir las ecuaciones del ecualizador IIR de

Kalman, cuya estructura se representa en la Figura 3:

ESTIMACIÓN

Ŝ(k/k) = Ŝ(k/k-1) + K(k)[x(k) – H(k) Ŝ(k/k-1)]

PREDICCIÓN

Ŝ(k/k-1) = aŜ(k-1/k-1)

GANANCIA

K(k) = V(k/k-1)HT(k)[H(k)V(k/k-1)HT(k) + n2(k)]-1

K(k) = [k0(k) k1(k) ... kd(k)]T

COVARIANZA DEL ERROR

V(k/k-1) = aV(k-1/k-1)aT(k) + bbTs2(k)

V(k/k) = [I – K(k)H(k)]V(k/k-1)]

La idea del montaje adaptativo (Figura 4) es que funcione de la siguiente forma. Habrá

un periodo inicial de entrenamiento en el cual se transmitirá una secuencia

predeterminada que servirá para que el algoritmo de identificación de sistemas

______________________________________________________________________

z-1 z-1z-1

0

h0 hM-1h1

k0 k1 kM-1

X(k)

Kd

Ŝ(k-d)

Figura 3:Estructura del Ecualizador (Tomada de [Mulgrew88])

________________________________________________________________TEAN

construya una estimación ĥ de la respuesta al impulso del canal h. A continuación se le

pasa esa estimación al Filtro de Kalman para que se inicialice, y comienza la

transmisión de datos. El hecho de utilizar una secuencia de entrenamiento supone que

la inicialización del Filtro será exacta.

El vector de estado en el instante cero estará constituido por las d últimas muestras de

la secuencia de entrenamiento; la matriz de covarianza de error será la matriz nula

puesto que el vector de estado se conoce con probabilidad uno.

El Filtro de Kalman proporciona la varianza del error de las estimaciones en la diagonal

principal de la matriz de covarianza de error. Una vez que el Filtro de Kalman ha

alcanzado su estado estacionario, se usará como entrada para el decisor el elemento

del vector de estado que alcance cierto umbral según el criterio que se haya decidido

utilizar.

Durante la transmisión de datos se puede utilizar el ecualizador adaptativo en la

modalidad de decisión directa (la otra posición de los conmutadores); posición en la

cual se usará la salida del decisor como entrada para el algoritmo de identificación de

______________________________________________________________________

Entrenamiento

Entrenamiento

X(k)s(k)FIR

CANAL

n(k)

m(k-d)ŝ(k-d)

DECISOR

ECUALIZADOR DE KALMAN

K

a

z--1

ĥ

—

FILTRADO ADAPTATIVO

LMS

IDENTIFICACIÓN DE SISTEMA

—

z--d

Figura 4:Montaje del Ecualizador Adaptativo (Tomada de [Mulgrew88])

________________________________________________________________TEAN

sistemas. Con esta forma de funcionamiento se supone que si el decisor toma las

decisiones correctas, el ecualizador será capaz de seguir las variaciones temporales del

canal. En esas condiciones la salida del decisor representará la entrada del canal con

un retardo de d muestras. De forma que para que el algoritmo de identificación de

sistemas funcione bien hay que meter un retardo de d muestras entre la salida del canal

y la entrada del algoritmo de identificación de sistemas.

Estos dibujos y esas ecuaciones son las que uno se encuentra en los libros de filtrado

adaptativo; y es lo que se mete en la máquina y funciona. La idea ahora es ir a los

artículos de R. E. Kalman donde se presentan los conceptos en los que se basan estas

estructuras y modelos...

Este hombre, R. E. Kalman, se propone en su artículo de 1960 abordar el problema del

filtrado y la predicción lineal usando la representación de Bode-Shannon de procesos

estocásticos y el método de transición de estados para el análisis del modelo de

______________________________________________________________________

________________________________________________________________TEAN

sistema dinámico que propone como solución al mismo. De hecho la aportación más

importante que hace este autor es la formulación del problema de Wiener desde un

punto de vista diferente al que aportaba la teoría clásica de análisis de sistemas: un

planteamiento en términos del espacio de estados; planteando también la dualidad

entre el problema de la estimación lineal óptima y el problema de control determinístico

del regulador sin ruido.

Según el planteamiento original que plantea el autor, se trata de obtener la

especificación de un sistema lineal dinámico para la predicción, separación o detección

de una señal aleatoria.

El artículo presenta un nuevo enfoque del problema de Wiener que ayuda a superar

algunas limitaciones de los métodos disponibles en la época que coartaban su

potencialidad y utilidad práctica en aplicaciones de ingeniería.

Hay una serie de suposiciones básicas de común acuerdo que están debajo de toda la

exposición de Kalman y el desarrollo de su modelo:

Todos los cálculos que se realizan y los resultados que se presentan se apoyan

en el conocimiento a priori de los estadísticos de primer y segundo orden.

Las señales aleatorias arbitrarias se representan como la salida de un sistema

lineal dinámico excitado por señales aleatorias incorreladas (la típica suposición

que se suele adoptar en el estudio de sistemas lineales que asume la existencia

de ruido blanco gausiano aditivo a la entrada del sistema provocando un

determinado error en las observaciones). Es en este sentido donde el autor ya

difiere del enfoque clásico para la descripción de los sistemas: Kalman propondrá

un sistema de ecuaciones en diferencias (o diferenciales en el caso contínuo) de

primer orden). El modelo se basa en los conceptos de ESTADO del sistema, y

TRANSICIÓN de estado.

La solución del problema de Wiener nos proporcionará los coeficientes del filtro

a partir de la matriz de covarianza.

Prácticamente la totalidad de su trabajo trata con sistemas discretos: que en los

términos en que los vamos a tratar se pueden considerar como sistemas

contínuos que han sido discretizados mediante muestreo.

______________________________________________________________________

________________________________________________________________TEAN

El autor formula sus resultados en forma de teoremas expresados en un lenguaje

matemático más o menos formal; utilizando una notación, que en mi opinión no es la

más afortunada que se le podía haber ocurrido.

De esos resultados que aparecen en los artículos, puede que aparezca alguno a lo

largo del documento, pero en la lectura que yo he hecho me ha parecido más

destacable el planteamiento que propone del problema de Wiener, los conceptos que

subyacen a todo el aparato estadístico que soporta sus razonamientos, y el alcance de

los mismos; no me gustaría que después de leer estos artículos sólo pudiera extraer un

par de algoritmos numéricos muy eficientes para meter en el MatLab, o cualquier otra

herramienta al uso, y jugar con ellos para sacar un montón de gráficas más o menos

vistosas donde se ponga de manifiesto, por ejemplo, la diferencia en la velocidad de

convergencia de los mismos.

Es por esto que nos fijaremos en las bases conceptuales que soportan el razonamiento

que sigue Kalman para construir su modelo. En este sentido, lo primero que podemos

destacar son una serie de consideraciones acerca del problema de la construcción de

estimaciones óptimas.

El problema parte de la suposición de que tenemos una señal x1(t) a la que se le suma

un ruido x2(t); partiremos por tanto de una serie temporal de observaciones y(t) = x1(t) +

x2(t), i.e. supondremos conocidos los valores de y(t0), ... , y(t). De lo que se trata es de

ver qué podemos inferir sobre los valores de la señal en un instante t1.

Este planteamiento nos lleva directamente a tratar con un problema de estadística y

teoría de la probabilidad, que se podría enunciar en los siguientes términos:

Dado un conjunto de observaciones de y(t) (una variable aleatoria V.A.), (t0), ... , (t)

deberíamos poder calcular la probabilidad de la ocurrencia simultánea de varios valores

(t) de la V.A. x1(t1), para lo cual propone echar mano de la función de probabilidad

condicional F(1) = Pr [x1(t1) 1|y(t0)=(t0), ... , y(t)=(t)].

La idea es que cualquier estimación será una función de esa distribución, es decir, una

función determinística de las V.A. y(t0), ... ,y(t) que podemos poner como X1(t1|t).

______________________________________________________________________

________________________________________________________________TEAN

Se definirá una función de pérdidas para penalizar las estimaciones incorrectas; y será

una función positiva y no decreciente del error cometido en la estimación = x1(t1) –

X1(t1). El criterio para elegir X1 será normalmente la minimización de las pérdidas

medias, o el “riesgo” en la estimación.

E{L[x1(t1) – X1(t1)]} = E[E{L[x1(t1) – X1(t1)]|y(t0), ... ,y(t)}] = E{L[x1(t1) – X1(t1)]|y(t0), ... ,y(t)}

Como concepto importante que soporta el razonamiento que sigue Kalman a la hora de

atacar el problema de la estimación lineal, aparece la proyección ortogonal.

La idea es que los resultados obtenidos mediante estimación lineal se pueden mejorar

con una estimación no lineal sólo cuando los procesos aleatorios son no gausianos y

además considerando funciones de distribución de tercer orden.

Para justificarse, Kalman busca una interpretación geométrica, y nos plantea las

siguientes consideraciones. Si tenemos las V.A. y(t0), ... ,y(t) el conjunto de todas las

combinaciones lineales (i=t0,...,t)aiy(i) formará un espacio vectorial (variedad lineal)

ℒ(t), que se puede considerar como un subespacio de dimensión finita del espacio

formado por todas las posibles observaciones.

En esas condiciones, dados dos vectores u,vℒ(t) diremos que son ortogonales sii Euv

= 0. Usando el método de Schmidt podremos seleccionar una base ortonormal de ℒ(t)

{et0, ... ,et} ; Eeiej = ij = 1 si i = j con i,j = t0, ... ,t

Χ = (i=t0,...,t)aiei, Χℒ(t)

Cualquier V.A. se podrá descomponer de forma unívoca en una parte Χ que caerá en

ℒ(t), y otra parte X’ que será ortogonal a ℒ(t)

X = Χ + X’ = (i=t0, ... ,t)(Exei)ei + X’

Con lo cual podremos escribir

E X’ei = E(x – Χ)ei = Exei - E Χ ei = 0 :(i=t0, ... ,t)

Puesto que las coordenadas son únicas.

Se dice entonces que Χ es la proyección ortogonal de x en ℒ(t).

A Χ se le puede ver como el vector de ℒ(t) (es decir una función lineal de las V.A.

y(t0), ... ,y(t)) que minimiza la función de péridas cuadrática. Eso quiere decir que si

cogemos cualquier otro vector de ℒ(t), w, tendremos

______________________________________________________________________

________________________________________________________________TEAN

E(x-w)2 = E(x’ + x – w)2 = E[(x –x) + (x – w)]2

Como resulta que x’ es ortogonal a todos los vectores de ℒ(t), y en particular a (x – w),

tendremos que

E(x-w)2 = E(x –x)2 + (x – w)2 E(x – x)2

Si además w minimiza la función cuadrática de pérdidas, E(x – w)2 = 0, y eso significa

que x = w

Lo que viene a decir este razonamiento es que si tenemos {x(t)} e {y(t)} procesos

aleatorios de media nula, y un conjunto de observaciones y(t0), ... ,y(t) podemos decir

que si ocurre que los procesos son gausianos, o bien que la estimación óptima se

restringe a una función lineal de las V.A. observadas y la función de pérdidas L() = 2

entonces x*(t1|t) la estimación óptima de x(t1) dados y(t0), ... ,y(t) coincide con x(t1|t) la

proyección ortogonal de x(t1) en ℒ(t). Esto lo expresa Kalman de una forma más

compacta como x(t1|t) = x*(t1|t) = Ê[x(t1)| ℒ(t)] ; que también se puede interpretar como

que en las condiciones mencionadas la mejor estimación que se puede hacer es una

combinación lineal de las observaciones previas. Estamos diciendo que el estimador

óptimo es un filtro al que le metemos los valores reales de las V.A. observables.

La aportación importante que hace Kalman a la resolución del problema del filtrado

lineal es su propuesta de un modelo para los procesos estocásticos que le servirá para

introducir un nuevo enfoque en el Análisis de Sistemas.

La aproximación habitual parte de la consideración de los procesos aleatorios

macroscópicos como procesos gausianos independientes.

Cuando tal aproximación no es aplicable se suele explicar la dependencia estadística

entre las observaciones de señales aleatorias en instantes temporales distintos con la

suposición de la existencia de un sistema dinámico entre la fuente y el observador.

Se puede ver por tanto una señal aleatoria como la salida de un sistema dinámico

excitado por una serie de procesos estocásticos gausianos independientes.

Con esto se trata de aprovechar la buena costumbre que tienen las señales aleatorias

gausianas de seguir siendo gausianas al pasar por un sistema lineal.

______________________________________________________________________

________________________________________________________________TEAN

Si hemos supuesto que las fuentes son gausianas e independientes, y las señales

aleatorias que forman el conjunto de observaciones son también gausianas, podremos

suponer que el sistema dinámico que hay entre la fuente y el observador es LINEAL.

Además sabemos que dado un proceso estocástico del que se conocen sus

estadísticos de primer y segundo orden, podemos encontrar un proceso gausiano con

las mismas propiedades.

Con estas ideas iniciales, Kalman se propone buscar la forma de describir el sistema en

base al concepto de ESTADO. El autor entiende por estado un conjunto de información

cuantitativa que constituye la cantidad mínima de datos que necesitamos conocer

acerca del comportamiento del sistema en el pasado para poder predecir su futuro.

La dinámica del sistema quedará descrita en términos de las TRANSICIONES entre

estados.

El sistema quedará descrito por un modelo como el siguiente:

dx/dt = F(t)x(t) + D(t)u(t)

y(t) = M(t)x(t)

donde x =[x1 ... xn]T es el estado del sistema, y además

u(t)ℳm1 ; mn

F(t)ℳnn ; D(t)ℳnm

y(t)ℳp1 ; M(t)ℳnp ; pn

Es un sistema en tiempo contínuo. Si suponemos que es un sistema estacionario (es

decir que todos los coeficientes de F, D y M son constantes: un sistema LTI), y que u(t)

es constante en el periodo de muestreo, u(t+) = u(t) ; 01 , t=0,1, ... podemos

transformarlo en un sistema discreto. Kalman lo expresa así:

X(t+1)=(1)x(t) + (1)u(t) ; t=0,1,...

(1) = exp F = (i=0..)(Fi/i!) ; F0=I (1) = (∫0,1expFd)D

Si no es estacionario, podremos ponerlo como

X(t+1)=(t+1;t)x(t) + (t)u(t) ; t=0,1,...

y(t) = M(t)x(t)

______________________________________________________________________

________________________________________________________________TEAN

(t+1;t) es la MATRIZ DE TRANSICIÓN, y es tal que:

(t2;t1) indica la transición desde t1 hacia t2

(t;t) = I

(t;s) (s;r) = (t;r)

-1(t;s) = (s;t)

Si el sistema es LTI, (t+1;t) = (t+1-t) = (1) = Constante y (t;) = exp F(t-)

Por otro lado, {u(t)} será un proceso vectorial, independiente, y gausiano de media nula:

Eu(t) = 0, t

Eu(t)u’(z) = 0 si tz

Eu(t)u’(t) = Q(t)

{x(t)} también será un proceso gausiano de media nula, pero no independiente tal que

x(t) = (r=- .. t-1) (t;r+1)u(r)

Si ts : Ex(t)x’(t) = (r=- .. s-1) (t;r+1)Q(r) ’(s;r+1)

Partiendo del modelo de sistema hemos llegado a la matriz de covarianzas Ex(t)x’(t); sin

embargo lo normal en los problemas del mundo real es que intentemos llegar a la

matriz de covarianzas partiendo de las observaciones, para luego construir el modelo a

partir de la matriz.

Con este modelo de sistema se puede ver el problema de Wiener de forma natural.

Considerando el siguiente sistema dinámico:

X(t+1) = (t+1;t)x(t) + u(t)

Y(t) = M(t)x(t)

Donde u(t) es un proceso gausiano independiente vectorial (vectores de n

componentes), y de media nula; x(t) es un vector de n componentes; y(t) tiene p

componentes (pn); (t+1;t)ℳnn ; M(t)ℳpn siendo sus elementos funciones

determinísticas del tiempo (señales).

Dadas las observaciones y(t0), ... ,y(t)

Encontrar un estimador x*(ti|t) de x(t1) que minimice la esperanza de la función de

pérdidas.

______________________________________________________________________

________________________________________________________________TEAN

Sabemos que la solución será la proyección ortogonal de x(t1) sobre la variedad lineal

generada por las V.A. observadas.

Al final, jugando con la variedad lineal generada por todas las observaciones y la

proyección ortogonal, Kalman propone la siguiente expresión para el Filtro Óptimo:

X*(t+1) = *(t+1;t)x*(t|t-1) + *(t)y(t)

*(t+1;t) = (t+1;t) - *(t)M(t)

cuyo diagrama de bloques podemos ver en la Figura 5.

En el artículo de Kalman, todo esto que caracterizado en forma de Teorema, donde se

recogen las relaciones de recursión que constituyen el algoritmo de resolución del

problema de Wiener:

La estimación óptima x*(t+1|t) de x(t+1) dados y(t0), ,,, ,y(t) se genera con el sistema

lineal dinámico

x*(t+1|t) = *(t+1;t)x*(t|t-1) + *(t)y(t)

El error de la estimación viene dado por

x(t+1|t) = *(t+1;t) x(t|t-1) + u(t)

______________________________________________________________________

(t+1;t)

(t+s;t+1)

*(t) M(t)

-I

x*(t|t-1)

x*(t+1|t-1)

x*(t+1|t)

x*(t+s|t)

y(t|t-1)ŷ(t|t-1)y(t)

hhfhhfhfhhhhhhhhhhhhhhhh

Figura 5:Diagrama de Bloques del Filtro Óptimo de Kalman

d

________________________________________________________________TEAN

La matriz de covarianza del error es

P*(t) = Ex(t|t-1)x(t|t-1)

La pérdida en media cuadrática es

(i=1, ... ,n)Exi2(t|t-1) = traza {P*(t)}

Las matrices *(t), *(t+1;t), y P*(t) siguen las siguientes relaciones de recursión

*(t) = (t+1;t)P*(t)M’(t)[M(t)P*(t)M’(t)]-1

*(t+1;t) = (t+1;t) - *(t)M(t)

P*(t+1) = *(t+1;t)P*(t)’(t+1;t) + Q(t)

Para obtener las iteraciones debemos especificar las covarianzas

P*(t0) de x(t0) y Q(t) de u(t)

La generalización para cualquier s0, si es inversible, vendrá dada por

x*(t+s|t) = (t+s;t+1)x*(t+1|t) = (t+s;t+1)*(t+1;t)(t;t+s-1)x*(t+s-1|t-1) + (t+s;t+1)*(t)y(t)

Si eliminamos * y * de las relaciones de recursión, obtendremos una ecuación en diferencias

no lineal

P*(t+1) = (t+1;t){P*(t) – P*(t)M’(t)[M(t)P*(t)M’(t)]-1P*(t)M(t)}’(t+1;t) + Q(t) ; tt0

Es interesante destacar en una primera lectura del artículo el planteamiento que hace el

autor del problema de control determinístico mostrando su dualidad con el de la

estimación lineal óptima.

Considerando el sistema dinámico

x(t+1) = (t+1;t)x(t) + M(t)u(t)

donde

x(t)ℳn1 ; u(t)ℳm1

(t)ℳnn ; M(t)ℳnm ; mn

siendo sus elementos funciones no aleatorias del tiempo

El problema que se plantea es el siguiente:

______________________________________________________________________

________________________________________________________________TEAN

Dado cualquier estado x(t) en un instante t, encontrar la secuencia u(t), ... ,u(T) de

vectores de control que minimicen el índice de rendimiento

V[x(t)] = (=t .. T+1)x’()Q()x()Donde Q(t) es una matriz definida positiva cuyos elementos son señales

determinísticas.

En este problema NO HAY CONSIDERACIONES ESTADÍSTICAS.

Kalman propone que la dualidad entre los problemas de estimación y control se puede

ver en el siguiente sentido:

Sea 0. Si en las relaciones del problema de estimación sustituimos cualquier matriz

de la forma X(t) = X(t0+) por X’(t) = X’(T-) obtendremos las relaciones

correspondientes al problema de control. Y a la inversa, si sustituimos cada matriz X(T-

) en las relaciones de control por X’(t0+), obtendremos las relaciones para la

estimación lineal óptima.

Matemáticamente lo que significa la dualidad es que ambos problemas se reducen a la

resolución de una ecuación no lineal del tipo de la de Wiener-Hopf, que es la que surgía

de la eliminación en las relaciones de recursión de dos de las matrices.

En la Figura 6 se puede ver una representación esquemática del problema del

regulador sin ruido tomada del artículo del propio Kalman

______________________________________________________________________

(t+1;t)

M(t) *(t)

-I

x(t+1)U*(t) hhfhhf

hfhhhhhhhhhhhhhhhh

Figura 6:Diagrama de Bloques del Regulador Sin Ruido

d

Sistema Físico a controlar

________________________________________________________________TEAN

Teniendo en cuenta que

U*(t) = - *(t)x(t)

Las ecuaciones para el control óptimo son

x(t+1) = * (t+1;t)x(t)

V*[x(t)] = x’(t)P*(t-1)x(t)

Las relaciones de recursión quedan así

* (t) = [M’(t)P*(t) M(t)]-1 M’(t) P*(t) (t+1;t)

* (t+1;t) = (t+1;t) - M(t) * (t)

P*(t+1) = ’ (t+1;t) P*(t) * (t+1;t) + Q(t)

Para tT

Siendo inicialmente P*(T) = Q(T+1)

Conclusiones

Los dos artículos de Kalman en los que se basa este documento son unos clásicos

dentro del vasto campo del filtrado y la predicción lineal; constituyeron la base para el

desarrollo de un serie de algoritmos de resolución del problema de Wiener como los

que se mencionaron más arriba y que son y han sido de amplia aplicación en los

campos de las comunicaciones y el control. Resulta clarificador estudiar el fundamento

conceptual que hay detrás de los algoritmos que metemos en las máquinas a la hora de

resolver los problemas reales; en este caso resulta interesante ver el jugo que se le ha

sacado a una forma intuitiva de representación del problema en base al espacio de

estados y la matriz de transición; el cómo tener un adecuado modelo de los procesos

estocásticos ayuda a llegar al corazón del problema de estimación dándole una

______________________________________________________________________

________________________________________________________________TEAN

interpretación geométrica a partir de la proyección ortogonal, que también da pie a

plantear la Ecuación Integral de Wiener-Hopf.

Sobre todo resulta estimulante pensar en la idea de la dualidad entre dos problemas

que a priori parecen no tener mucho que ver: un problema, el de la ESTIMACIÓN (ya

sea interpolación, filtrado o predicción), cubierto por un aparato estadístico que puede

conseguir desviar nuestra atención del fondo de la cuestión; y otro, el del CONTROL

determinístico de una planta sin ruido donde no se hace ninguna consideración

probabilística. Una dualidad que matemáticamente nos lleva a la ecuación de Wiener-

Hopf, y que no encontraba significado físico bajo el enfoque de Kalman, que había

estudiado ambos problemas por separado.

Bibliografía

Para la elaboración de este documento se han consultado, aparte de los apuntes de la

asignatura, las siguientes fuentes

1. Mulgrew, Bernard and F.N. Cowan, Colin “Adaptive Filters and Equalisers” Kluwer

Academic Publishers USA, 1988

2. Wiener, Norbert “CIBERNÉTICA o el control y comunicación en animales y

máquinas” 2ª edición Tusquets Editores, Barcelona 1998

3. R. E. Kalman “A New Approach to Linear Filtering and Prediction Problems” Trans.

ASME, Journal of Basic Engineering pp. 35-45, (March 1960)

______________________________________________________________________

________________________________________________________________TEAN

4. R. E. Kalman and R. S. Bucy “ New Results in Linear Filtering and Prediction Theory”

Trans. ASME, Journal of Basic Engineering pp. 95-108, (March 1961)

5. R. E. Lawrence and Kaufman “The Kalman Filter for the Equalization of a Digital

Communications Channel” IEEE Transactions on Communications Technology vol.

19(6) pp. 1137-1141, (Dec. 1971)

______________________________________________________________________

________________________________________________________________TEAN

Temas Especiales de Análisis Numérico

José Manuel Corrales Sendino

Curso 1999/2000

Primer Semestre

______________________________________________________________________