Introducción (planteamiento)
Dado el siguiente modelo lineal, estimar el único estado (constante) utilizando un filtro de Kalman a partir de la medición de una señal tipo escalera con ruido (con segmentos constantes cuyos valores se deben estimar) y analizar los resultados:
{˙x=0+wy=x+v
Método (plan de solución)
Pasos: (i) implementar en MATLAB el filtro de Kalman; (ii) implementar una señal de salida tipo escalera y adicionar un ruido blanco gaussiano y un ruido coloreado; (iii) estimar el valor del estado utilizando diferentes valores de [P(0),Q,R]; (iv) analizar los resultados.
Resultados (solución)
Modelo discretizado utilizando el método de Euler:
{x(k+1)=x(k)y(k)=x(k)
Filtro de Kalman con G=I y matrices equivalentes a escalares por ser un sistema de primer orden:
(1) K(k)=P(k|k−1)[P(k|k−1)+R]−1(2) ˆx(k|k)=ˆx(k|k−1)+K(k)[y(k)−ˆx(k|k−1)](3) P(k|k)=[I−K(k)]P(k|k−1)(4) ˆx(k|k−1)=ˆx(k−1|k−1)(5) P(k|k−1)=P(k−1|k−1)+Q
Código de MATLAB que implementa el filtro de Kalman para el problema propuesto, adicionando un ruido coloreado a la salida:
Ts = 0.01;N = 1000;t = (0:Ts:(N-1)*Ts)';rng(1)e = randn(N,1);H = tf([1 0.1],[1 0.5],'Ts',Ts);v = lsim(H,e,t);Av = 0.1; % Desviación estándar de la perturbacióny1 = [0.4267*ones(1,N/4) 0.9372*ones(1,N/4) -0.1298*ones(1,N/4) -0.7833*ones(1,N/4)]';y = y1 + (Av/std(v))*v; % Señal con ruido coloreadox0 = 0;P0 = 1;Q = 0.001;R = Av^2; % Condiciones iniciales y parámetros de diseñox_act = zeros(N,1);P_act = zeros(N,1);K = zeros(N,1);x_pre = x0;P_pre = P0;for k=1:N % Filtro de Kalman% Corrección del estadoK(k) = P_pre/(P_pre + R);x_act(k) = x_pre + K(k)*( y(k) - x_pre );P_act(k) = ( 1 - K(k) )*P_pre;% Actualización del estadox_pre = x_act(k);P_pre = P_act(k) + Q;ende_est = y-x_act; % Residuosfigure(1)sgt = sgtitle(['x(0) = ' num2str(x0) ', P(0) = ' num2str(P0) ', Q = ' num2str(Q) ', R = ' num2str(R)]);sgt.FontSize = 10;subplot(2,2,1) % Estado estimado con la desviación estándar a partir de P_actplot(t,y,'c',t,y1,'r',t,x_act+sqrt(P_act),'r--',t,x_act-sqrt(P_act),'r--',t,x_act,'k'), xlabel('t (seg)'), title('x(t)')subplot(2,2,2)plot(t,e_est,t,zeros(N,1)), xlabel('t (seg)'), title('e(t)')subplot(2,2,3)plot(t,P_act), xlabel('t (seg)'), title('P(t)')subplot(2,2,4)plot(t,K), xlabel('t (seg)'), title('K(t)')figure(2)whiteness_test(e_est)
Las siguientes figuras muestran los resultados con x(0) = 0 y diversos valores de {P(0), Q, R} (ver el título de cada figura).
Figura 1 con una confianza errónea en el estado inicial debido a P(0) = 0 y con una mayor confianza en la predicción que en la medición (Q < R), donde se observa un error inicial en la estimación, una rápida convergencia de la estimación con cada cambio de la salida medida y un valor de K menor que 0.3:
Figura 2 con una menor confianza en la estimación que en la medición (Q > R), donde se observa una rápida convergencia de la estimación con cada cambio de la salida medida, una varianza grande del estado estimado y un valor de K cercano a 1:
Figura 3 con poca confianza en el estado inicial debido a P(0) = 1 y con una mayor confianza en la predicción que en la medición (Q < R), aunque con un valor de R mayor al real, donde no se observa un error de estimación al inicio y hay una menor varianza de la estimación con una mayor incertidumbre:
Figura 4 con una mayor confianza en la predicción que en la medición (Q < R), aunque con un valor de R mucho mayor al real, donde se observa una menor varianza de la estimación, pero con una incertidumbre mucho mayor que el caso con R = 0.1, y con una baja velocidad de convergencia:
Figura 5 con una destacable mayor confianza en la predicción que en la medición (Q << R), donde se observa una también destacable menor varianza de la estimación con una incertidumbre mucho menor que en el caso con Q = 0.001 y una casi nula corrección (valor de K muy pequeño), pero con una baja velocidad de convergencia:
Figura 6 con la prueba de blancura con [P(0), Q, R] = [1, 0.00001, 0.01]:
Figura 7 con la prueba de blancura con [P(0), Q, R] = [1, 0.0001, 0.01], adicionando un ruido blanco gaussiano a la salida (v = e) con una mayor desviación estándar de 0.5, donde se observan mejores resultados en la prueba de blancura (aunque no significativos):
Discusión y verificación
Este ejercicio con matrices escalares muestra de manera simple la manera como se implementa un filtro de Kalman y cómo se comportan las matrices P y K durante la estimación, aunque no se deben generalizar los resultados a cualquier problema. Las pruebas se realizan con un ruido coloreado adicionado a la salida medida, sin errores significativos en las estimaciones; solo la figura 7 con un ruido blanco gaussiano muestra una mejor prueba de blancura.
Si no se tiene confianza en el estado inicial, lo recomendable es iniciar con un valor grande de la matriz P (comparar la figura 1 con las demás figuras).
Parece una buena idea utilizar valores de Q menores que R para reducir la varianza (comparar las figuras 1, 3, 4 y 5 con la figura 2), aunque valores de Q mucho más pequeños que de R hacen más lenta la convergencia del algoritmo (figuras 4 y 5).
Hacer Q mayor que P hace que la varianza aumente considerablemente en la estimación.
Dejar R en su valor experimental y ajustar Q parece ser una buena idea: la figura 1 es bastante aceptable corrigiendo el estado inicial, con poca varianza y buena velocidad de convergencia; si se sigue disminuyendo Q se obtiene una baja velocidad de convergencia (figura 5). En este caso los residuos también tienen un buen comportamiento (figura 6) si no se tienen en cuenta los cambios bruscos en la salida medida, lo cual no es habitual en los procesos.
Finalmente, se observa en todos los casos la convergencia de P y K a un valor, lo cual demuestra que el ruido corresponde a un proceso estocástico estacionario. De esta manera, tener una información de P y K ayuda a predecir de manera rápida el estado sin necesidad de partir siempre de un desconocimiento total del comportamiento pasado del sistema.
Comentarios