Práctica 4: Sensores Inteligentes

Observador de Estados

1. Aprender a utilizar herramientas matemáticas (Matriz de Observabilidad y Rango) para determinar si es teóricamente posible reconstruir los estados internos de un sistema a partir de sus salidas medidas.
2. Construir un estimador de estados asintótico para un sistema dinámico (circuito RLC), calculando la matriz de ganancia L mediante la asignación de polos para asegurar la convergencia del error a cero.
3. Evaluar experimentalmente cómo la ubicación de los polos del observador afecta tanto la velocidad de convergencia de la estimación como la sensibilidad al ruido de medición.
4. Contrastar el diseño manual de polos (mediante el comando place) frente a un diseño óptimo considerando el ruido estocástico (Observador LQR o Filtro de Kalman Estacionario), para entender las ventajas de este último en entornos ruidosos.

Marco Teórico: Observadores de Estado.

Definición y Propósito

En el control de sistemas dinámicos, a menudo es necesario conocer el vector de estado completo para implementar leyes de control (como la retroalimentación de estados). Sin embargo, frecuentemente no es posible medir todos los estados directamente debido a limitaciones físicas o costos de sensores.

Un observador de estados es un sistema dinámico construido matemáticamente que estima el vector de estados basándose en las variables que sí son medibles: la entrada de la planta y la salida de la planta .

Donde la entrada de un sensor entra a la estructura matemática descrita en la Figura 1, podemos estimar el estado . Por otro lado, si definimos el error de estimación como , la dinámica del error resulta en

Para solucionar las limitaciones del lazo abierto, se introduce un término de retroalimentación que corrige la estimación comparando la salida real con la salida estimada , como se muestra en la Figura 2.

Donde, podemos ver de forma directa, que la dinámica de convergencia de depende de los valores característicos de la matriz . Por lo que como se indica en [1], si el par es observable, podemos elegir una matriz Ltal que los valores propios estén en cualquier ubicación deseada, lo que permite:

1. Asegurar que el error converja a cero.
2. Hacer que el error converja mucho más rápido que la dinámica de la planta.

Observabilidad

Un sistema es observable si el estado inicial puede ser determinado a partir de la historia de la entrada y la salida en un intervalo de tiempo finito [3]. Lo anterior puede entenderse como si una variable de estado no tiene ningún efecto (directo o indirecto) sobre la salida, no puede ser "observada" ni estimada.

La prueba para saber si un sistema es observables, es mediante la Matriz de Observabilidad , la cual se define como:

Paso 1. Definición del sistema y prueba de Observabilidad.

Dado que el sistema es completamente observable, pasemos con el diseño de la ganacia del Observador

• ¿Qué realiza el comando 'obsv'?
• Investigue el teorema de dualidad.

Paso 2. Diseño de la Ganancia del Observador.

Recordemos que la dinámica del error está gobernada por la matriz , como se indica en la ecuación (7). Elegir una L tal que los valores propios de sean estables y rápidos.

• Investigue la restricción del comando 'place' al indicar los polos deseados.
• Investigue y modifique el código anterior para utilizar el comando 'acker'.

t_span = [0 4];

 

% Definir Entrada u(t)

u_func = @(t) 12 * (sin(2*pi*0.4*t) > 0);

 

% Pre-generar el Ruido

sigma = 0;% Magnitud del ruido. Darle un valor >0 para simular con ruido <1

fs_ruido = 2000; % Alta frecuencia para que parezca continuo

t_ruido = -0.01 : (1/fs_ruido) : (t_span(end)+0.01);

ruido_vals = sigma * randn(size(t_ruido));

% Interpolador para el ruido

ruido_func = @(t) interp1(t_ruido, ruido_vals, t, 'previous');

 

% Condiciones Iniciales

x0_real = [0; 0];

x0_est = [0.15; 10];

z0 = [x0_real; x0_est];

 

% Ejecutar ODE45 con la función de ruido pre-generado

[t, z] = ode45(@(t,z) dynamics_rlc(t, z, A, B, C_mat, L_obs, u_func, ruido_func), t_span, z0);

% Desempaquetar resultados

x_real = z(:, 1:2);

x_est = z(:, 3:4);

 

% Recalcular vectores exactos para graficar

u_vec = arrayfun(u_func, t); % Entrada exacta

ruido_vec = arrayfun(ruido_func, t); % Ruido exacto

y_clean = (C_mat * x_real')'; % Salida limpia

y_medida = y_clean + ruido_vec; % Salida ruidosa (Lo que entró al observador)

 

 

% ENTRADA DEL SISTEMA

figure('Name', 'Entrada del Sistema', 'Color', 'w');

plot(t, u_vec, 'b', 'LineWidth', 2);

ax = gca;

ax.Color = 'w';

ax.XColor = 'k';

ax.YColor = 'k';

ax.GridColor = 'k';

ax.GridAlpha = 0.15;

title('Entrada del Sistema - Voltaje (u(t))', 'FontSize', 10, 'FontWeight', 'bold', 'Color', 'k');

ylabel('Voltaje (V)', 'FontSize', 10, 'Color', 'k');

xlabel('Tiempo (s)', 'FontSize', 10, 'Color', 'k');

grid on;

xlim([0 t_span(end)]);

ylim([min(u_vec)-2 max(u_vec)+2])

lgd = legend('Voltaje de Entrada', 'Location', 'best');

set(lgd, 'TextColor', 'k', 'Color', 'w', 'EdgeColor', 'k', 'FontSize', 10);

 

% SALIDA DEL SISTEMA

figure('Name', 'Salida del Sistema', 'Color', 'w', 'Position', [100 100 800 400]);

plot(t, y_medida, 'k-', 'LineWidth', 0.5, 'Color', [0 0 0 0.3]); hold on;

plot(t, y_clean, 'r', 'LineWidth', 2);

% Ruido en gris tenue

ax = gca;

ax.Color = 'w';

ax.XColor = 'k';

ax.YColor = 'k';

ax.GridColor = 'k';

ax.GridAlpha = 0.15;

title('Salida del Sistema - Medición (y(t))', 'FontSize', 10, 'FontWeight', 'bold', 'Color', 'k');

ylabel('Salida', 'FontSize', 10, 'Color', 'k');

xlabel('Tiempo (s)', 'FontSize', 10, 'Color', 'k');

grid on;

xlim([0 t_span(end)]);

lgd = legend('Medición', 'Salida Real', 'Location', 'best');

set(lgd, 'TextColor', 'k', 'Color', 'w', 'EdgeColor', 'k', 'FontSize', 10);

 

 

% 3. CORRIENTE DEL INDUCTOR

figure('Name', 'Corriente Inductor', 'Color', 'w', 'Position', [100 100 800 500]);

plot(t, x_est(:,1), 'b--', 'LineWidth', 2); hold on;

plot(t, x_real(:,1), 'r', 'LineWidth', 2);

ax = gca;

ax.Color = 'w';

ax.XColor = 'k';

ax.YColor = 'k';

ax.GridColor = 'k';

ax.GridAlpha = 0.15;

title('Estimación de Corriente Inductor', 'FontSize', 10, 'FontWeight', 'bold', 'Color', 'k');

ylabel('Amperes (A)', 'FontSize', 10, 'Color', 'k');

xlabel('Tiempo (s)', 'FontSize', 10, 'Color', 'k');

lgd = legend('Corriente Estimada (Observador)','Corriente Real', 'Location', 'best');

set(lgd, 'TextColor', 'k', 'Color', 'w', 'EdgeColor', 'k', 'FontSize', 10);

grid on;

xlim([0 t_span(end)]);

 

 

% 4. VOLTAJE DEL CAPACITOR

figure('Name', 'Voltaje Capacitor', 'Color', 'w', 'Position', [100 100 800 500]);

plot(t, x_est(:,2), 'b--', 'LineWidth', 2); hold on;

plot(t, x_real(:,2), 'r', 'LineWidth', 2);

ax = gca;

ax.Color = 'w';

ax.XColor = 'k';

ax.YColor = 'k';

ax.GridColor = 'k';

ax.GridAlpha = 0.15;

title('Estimación de Voltaje Capacitor', 'FontSize', 10, 'FontWeight', 'bold', 'Color', 'k');

ylabel('Voltios (V)', 'FontSize', 10, 'Color', 'k');

xlabel('Tiempo (s)', 'FontSize', 10, 'Color', 'k');

lgd = legend('Voltaje Estimado (Observador)','Voltaje Real', 'Location', 'best');

set(lgd, 'TextColor', 'k', 'Color', 'w', 'EdgeColor', 'k', 'FontSize', 10);

grid on;

xlim([0 t_span(end)]);

 

 

% 5. ERROR ABSOLUTO - CORRIENTE

figure('Name', 'Error Absoluto - Corriente', 'Color', 'w', 'Position', [100 100 800 400]);

error_iL = abs(x_real(:,1) - x_est(:,1));

area(t, error_iL, 'FaceColor', [0.85 0.33 0.1], 'FaceAlpha', 0.4, 'EdgeColor', 'r', 'LineWidth', 1);

ax = gca;

ax.Color = 'w';

ax.XColor = 'k';

ax.YColor = 'k';

ax.GridColor = 'k';

ax.GridAlpha = 0.15;

title('Error Absoluto de Estimación - Corriente Inductor', 'FontSize', 10, 'FontWeight', 'bold', 'Color', 'k');

ylabel('|Error| (A)', 'FontSize', 10, 'Color', 'k');

xlabel('Tiempo (s)', 'FontSize', 10, 'Color', 'k');

grid on;

xlim([0 t_span(end)]);

error_max_iL = max(error_iL);

error_final_iL = error_iL(end);

text(t_span(end)*0.4, error_max_iL*0.6, ...

sprintf('Error máximo: %.4f A\nError final: %.4f A', error_max_iL, error_final_iL), ...

'FontSize', 10, 'Color', 'k', 'FontWeight', 'bold', 'BackgroundColor', 'w');

 

 

% 6. ERROR ABSOLUTO - VOLTAJE

figure('Name', 'Error Absoluto - Voltaje', 'Color', 'w', 'Position', [100 100 800 400]);

error_vC = abs(x_real(:,2) - x_est(:,2));

area(t, error_vC, 'FaceColor', [0.85 0.33 0.1], 'FaceAlpha', 0.4, 'EdgeColor', 'r', 'LineWidth', 1);

ax = gca;

ax.Color = 'w';

ax.XColor = 'k';

ax.YColor = 'k';

ax.GridColor = 'k';

ax.GridAlpha = 0.15;

title('Error Absoluto de Estimación - Voltaje Capacitor', 'FontSize', 10, 'FontWeight', 'bold', 'Color', 'k');

ylabel('|Error| (V)', 'FontSize', 10, 'Color', 'k');

xlabel('Tiempo (s)', 'FontSize', 10, 'Color', 'k');

grid on;

xlim([0 t_span(end)]);

error_max_vC = max(error_vC);

error_final_vC = error_vC(end);

text(t_span(end)*0.4, error_max_vC*0.6, ...

sprintf('Error máximo: %.4f V\nError final: %.4f V', error_max_vC, error_final_vC), ...

'FontSize', 10, 'Color', 'k', 'FontWeight', 'bold', 'BackgroundColor', 'w');

 

 

%VISUALIZACIÓN DEL RUIDO

figure('Name', 'Visualización del Ruido', 'Color', 'w', 'Position', [100 100 800 300]);

plot(t, ruido_vec, 'Color', [0.4 0.4 0.4], 'LineWidth', 1);

ax = gca;

ax.Color = 'w';

ax.XColor = 'k';

ax.YColor = 'k';

ax.GridColor = 'k';

ax.GridAlpha = 0.15;

title('Ruido Inyectado en el Sensor', 'FontSize', 10, 'FontWeight', 'bold', 'Color', 'k');

ylabel('Amplitud Ruido', 'FontSize', 10, 'Color', 'k');

xlabel('Tiempo (s)', 'FontSize', 10, 'Color', 'k');

grid on;

xlim([0 t_span(end)]);

 

 

%% FUNCIÓN DINÁMICA

function dz = dynamics_rlc(t, z, A, B, C, L, u_func, noise_func)

x = z(1:2);

x_hat = z(3:4);

 

u = u_func(t); % Evaluar entrada

ruido = noise_func(t); % Evaluar ruido pre-generado

 

% Sistema Real

dx = A * x + B * u;

 

% Medición ruidosa

y_medida = C * x + ruido;

 

% Observador

y_hat = C * x_hat;

dx_hat = A * x_hat + B * u + L * (y_medida - y_hat);

 

dz = [dx; dx_hat];

end

 

Actividad

1. Cambia N_obs de 2 a 0.5, 1, 5, y 10. Anota tus observaciones.
2. Cambia la variable sigma de 0 a 0.2 (o mayor) y N_obs en 2. Ejecuta la simulación. Observa la gráfica "Salida del Sistema". ¿Qué tan difícil es ahora distinguir la señal cuadrada original dentro de la mancha gris de ruido? ¿Crees que el observador pueda seguir funcionando bien?.
3. Mantén sigma = 0.2 Cambia N_obs de 2 a 10 (haciendo los polos 10 veces más rápidos/lejanos). Mira la línea azul punteada en "Corriente Inductor". ¿La estimación se volvió más limpia o más sucia? ¿Por qué al pedir rapidez ganamos ruido?
4. Con sigma = 0.2, cambia ahora N_obs a 0.5 (haciendo al observador más lento que el sistema físico). Observa la gráfica de "Error Absoluto". ¿Desapareció el ruido de alta frecuencia? Pero, ¿qué pasó con el "pico" de error al inicio (transitorio)? ¿Tarda más o menos en converger?
5. Investigue qué es un Observador LQR, comúnmente conocido en la literatura como Filtro de Kalman Estacionario.
6. ¿Cuál es la diferencia conceptual entre asignar los polos manualmente (como hicimos con place) y usar un criterio LQR?. ¿Qué representan las matrices de covarianza Q (ruido de proceso) y R (ruido de medición) en este contexto?
7. En su código, comente la línea donde usa place (% L_obs = ...) y copie el siguiente bloque de código para implementar el Filtro de Kalman Continuo.

Referencias