Práctica 7: Sensores Inteligentes
Sensor Virtual Térmico basado en Filtro de Kalman
Objetivos:
- Representar un sistema de parámetros distribuidos (ecuación de calor en una barra) utilizando el método de diferencias finitas para obtener un modelo en espacio de estados.
- Generar datos sintéticos que representen la temperatura real de la barra y añadir ruido gaussiano para simular un sensor físico real.
- Implementar un Filtro de Kalman discreto capaz de estimar la temperatura en puntos donde no existen sensores físicos.
- Evaluar cómo afectan las matrices de covarianza (Q y R) y las condiciones iniciales en la velocidad y precisión de la estimación.
Motivación:
En la industria moderna y la investigación (como en el monitoreo de baterías de litio, hornos industriales o procesadores de alto rendimiento), a menudo es imposible colocar sensores físicos en todas partes debido a limitaciones de espacio, temperaturas extremas o costos elevados.
Esto crea un problema: necesitamos conocer qué pasa dentro del sistema, pero solo podemos medir lo que pasa en la frontera. La motivación de esta práctica es demostrar cómo podemos usar el software (algoritmos matemáticos) para reemplazar al hardware. Al combinar un modelo físico con una sola medición barata en el exterior, podemos construir un "Sensor Virtual" que nos permite "ver" el interior del material como si tuviéramos cientos de sensores instalados.
Procedimiento:
A continuación, irás construyendo el sistema paso a paso. Lee la explicación y transcribe el código correspondiente en tu archivo de trabajo
Marco Teórico
Instrumentación en Sistemas Distribuidos.
La conducción de calor es un proceso fundamental de transferencia de energía térmica en medios sólidos, líquidos y gaseosos, gobernada por gradientes de temperatura, tal como se explica en la literatura clásica de [1]. Este fenómeno es crítico en aplicaciones que van desde sistemas industriales (hornos, intercambiadores) hasta dispositivos biomédicos y electrónicos [2].
En entornos experimentales, conocer el gradiente de temperatura suele requerir múltiples sensores distribuidos. Sin embargo, esta aproximación tiene limitaciones: los sensores pueden ser costosos, invasivos o técnicamente inviables en geometrías complejas. Esto da lugar a los llamados problemas inversos dinámicos, donde es necesario reconstruir información interna a partir de observaciones parciales y ruidosas en la frontera [2].
Sistemas de Parámetros Distribuidos (DPS).
En la ingeniería de control clásica, es común trabajar con sistemas de parámetros concentrados, donde se asume que las variables dependen únicamente del tiempo y se modelan mediante Ecuaciones Diferenciales Ordinarias (EDO). Sin embargo, fenómenos físicos como la transferencia de calor o el flujo de fluidos presentan variaciones tanto temporales como espaciales. Estos sistemas se denominan Sistemas de Parámetros Distribuidos (DPS) y se describen matemáticamente mediante Ecuaciones en Derivadas Parciales (EDP) [3].
Un ejemplo canónico es la ecuación de difusión de calor unidimensional para una barra homogénea de longitud L:
Donde α es la difusividad térmica. A diferencia de un sistema mecánico simple, aquí el "estado" es una función definida en todo el dominio espacial. Para este experimento, asumimos que los extremos de la barra se mantienen a temperatura ambiente (0°C), lo que implica que el calor se disipa por los bordes hacia el exterior, provocando que la temperatura del sistema decaiga eventualmente a cero.
Paso 1. Discretización Espacial (Método de Diferencias Finitas)
Como la computadora no puede resolver para infinitos puntos, dividimos la barra en Nnodos separados por una distancia 
. La segunda derivada espacial se aproxima numéricamente usando los vecinos del nodo actual [4]:
. La segunda derivada espacial se aproxima numéricamente usando los vecinos del nodo actual [4]:
Esto nos permite reescribir la ecuación parcial como un sistema lineal de ecuaciones diferenciales ordinarias (Espacio de Estados Continuo):
Donde la matriz Atiene una estructura tridiagonal característica 
escalada por 
, y B es la matriz de entrada, donde 
es el ruido de proceso térmico.
escalada por , y B es la matriz de entrada, donde es el ruido de proceso térmico.Transcriba la siguiente sección de código en su archivo .xml
clear all; close all; clc;
% Parámetros Físicos
L = 1; % Longitud de la barra (metros)
alpha = 0.01; % Difusividad térmica (m^2/s)
N = 20; % Número de nodos internos - Orden del modelo
% Parámetros de Discretización
dx = L / (N+1); % Distancia entre nodos (Delta x)
k_diff = alpha / dx^2; % Coeficiente de difusión numérico
% Construcción del modelo en espacio de estados
e = ones(N, 1);
A = spdiags([e -2*e e], -1:1, N, N);
A = k_diff * full(A); % Escalamos por la física
B = eye(N);
C = zeros(1, N);
C(1) = 1; % Solo tenemos sensor en el NODO 1
D = 0;
% Creación del objeto de sistema continuo
disp('Modelo continuo');
sys_cont = ss(A, B, C, D)
disp('Eigenvalores:');
disp(eig(A));
Preguntas de la sección:
- Observe los eigenvalores. ¿Son todos negativos? ¿Qué pasaria físicamente si hubiera un eigenvalor positivo?
- Si aumentamos
de 20 a 100, ¿qué le sucede al valor de 
y, por consecuencia, a los elementos de la matriz 
?¿La simulación se vuelve más lenta o más rápida?
Paso 2. Discretización Temporal.
Los sensores digitales no leen continuamente; leen cada cierto tiempo 
(tiempo de muestreo). Debemos convertir el sistema continuo a discreto:
(tiempo de muestreo). Debemos convertir el sistema continuo a discreto:
Donde Fes la matriz de transición de estado. Matemáticamente, la solución exacta es la exponencial matricial
Transcriba la siguiente sección de código en su archivo .xml
dt = 0.05; % Tiempo de muestreo (50 ms)
% Discretización "exacta" por Retenedor de orden zero
sys_d = c2d(sys_cont, dt, 'zoh')
F = sys_d.A; % Matriz de transición de estado discretizada
H = sys_d.C; % Matriz de medición discretizada
G = sys_d.B; % Matriz de entrada de ruido
Preguntas de la sección:
- Investigue los parámetros del comando c2d y que métodos de discretización puede realizar.
- Si el fenómeno de difusión es muy rápido (
muy grande) y su tiempo de muestreo grande (10 segundos por ejemplo), ¿podrá el sensor capturar los cambios? - Cambie el tiempo de muestreo a 100 segundos. Vuelva a calcular la matriz de transición de estado discreta y responda. ¿La matriz se parece a la identidad o cambia drásticamente? ¿Qué implica esto para la predicción del futuro?
Paso 3. Simulación.
Para validar el sensor virtual, necesitamos conocer los valores reales de temperatura contra la cual comparar nuestras estimaciones. En una aplicación real, esta verdad sería desconocida (por eso necesitamos el estimador), pero en simulación podemos generarla.
En lugar de un calentamiento uniforme, simularemos un escenario dinámico más complejo: un pico de calor concentrado en el centro de la barra [2]. Matemáticamente, esto se modela mediante una distribución Gaussiana como condición inicial:
Donde z es la posición a lo largo de la barra, 
la posición donde esta centrada la distribución, 
es la temperatura máxima y β define qué tan concentrado está el calor. Físicamente, esto equivale a haber calentado el centro de la barra con una fuente puntual justo antes de iniciar la prueba (
).
la posición donde esta centrada la distribución, es la temperatura máxima y β define qué tan concentrado está el calor. Físicamente, esto equivale a haber calentado el centro de la barra con una fuente puntual justo antes de iniciar la prueba ().En la práctica, ningún sensor es perfecto. Las lecturas de los termopares o RTDs siempre están contaminadas por ruido eléctrico, interferencias electromagnéticas o errores de cuantización. Para probar la robustez del Filtro de Kalman, modelamos esta imperfección añadiendo un término estocástico a la medición ideal
donde 
es la medición ruidosa que recibirá el microcontrolador, 
es la temperatura real en el punto de medición y 
es el ruido de medición, modelado como ruido blanco Gaussiano con desviación estándar 
es la medición ruidosa que recibirá el microcontrolador, es la temperatura real en el punto de medición y es el ruido de medición, modelado como ruido blanco Gaussiano con desviación estándar T_total = 20; % Tiempo total de simulación
t = 0:dt:T_total; % Vector de tiempo
steps = length(t);
% Temperatura maxima en C
T_inicial = 30;
%x0 = T_inicial * ones(N, 1);
x_pos = linspace(0, L, N);
x0 = 30 *exp(-10 * (x_pos - L/2).^2)';
% Simulación de la dinámica
[y_clean, t, x_true] = lsim(sys_d, zeros(steps, N), t, x0);
x_true = x_true';
% Añadir Ruido al Sensor Físico
sigma_sensor = 1.5; % Desviación estándar del ruido (grados Celsius)
ruido = sigma_sensor * randn(1, steps);
y_medida = y_clean' + ruido;
% VISUALIZACIÓN 2D
figure('Color', 'w');
imagesc(t, x_pos, x_true);
set(gca, 'YDir', 'normal');
colormap('jet');
c = colorbar;
c.Label.String = 'Temperatura (°C)';
set(c, 'Color', 'k');
set(c.Label, 'Color', 'k');
set(gca, 'FontSize', 10, 'LineWidth', 1.5, ...
'XColor', 'k', ...
'YColor', 'k');
xlabel('Tiempo (s)', 'FontSize', 12, 'Color', 'k');
ylabel('Posición en la Barra (m)', 'FontSize', 12, 'Color', 'k');
title('Evolución Espacio-Temporal', 'FontSize', 14, 'Color', 'k');
% VISUALIZACIÓN 3D
figure('Color', 'w', 'Position', [100 100 1000 600]);
[T_grid, X_grid] = meshgrid(t, x_pos);
surf(T_grid, X_grid, x_true);
shading interp;
colormap jet;
cb = colorbar;
set(cb, 'Color', 'k');
xlabel('Tiempo (s)', 'FontSize', 12, 'FontWeight', 'bold', 'Color', 'k');
ylabel('Posición (m)', 'FontSize', 12, 'FontWeight', 'bold', 'Color', 'k');
zlabel('Temperatura (°C)', 'FontSize', 12, 'FontWeight', 'bold', 'Color', 'k');
axis tight;
view(-45, 30);
grid on;
set(gca, 'Color', 'w', ...
'XColor', 'k', ...
'YColor', 'k', ...
'ZColor', 'k', ...
'GridColor', 'k', ...
'GridAlpha', 0.2);
% LECTURAS DEL SENSOR
figure('Color', 'w');
plot(t, x_true(1,:), 'k', 'LineWidth', 1.5); hold on;
plot(t, y_medida, 'r.', 'LineWidth', 2)
set(gca, 'Color', 'w', ...
'XColor', 'k', ...
'YColor', 'k', ...
'GridColor', 'k');
lgd = legend('Temperatura real', 'Medición ruidosa');
xlabel('Tiempo (s)')
ylabel('Temperatura (°C)')
set(lgd, 'TextColor', 'k', 'Color', 'w', 'EdgeColor', 'k');
Preguntas de la sección:
- Modifique sigma_sensor a 5. ¿Distingue la señal detrás del ruido? ¿Cree que el filtro a implementar podrá funcionar con tanto ruido?
- La condición inicial x0 pone calor en el centro. Sin embargo, nuestro sensor esta en x=0. Observe la gráfica y responda: ¿hay un retardo entre t=0 y cuando la temperatura empieza a subir en el sensor?¿Por qué?
Paso 4. Sensor virtual.
En la instrumentación industrial clásica, un sensor RTD (Resistance Temperature Detector) ubicado en el extremo de la barra (
) solo nos daría información puntual de esa posición y sería un error asumir que toda la barra se encuentra a la misma temperatura. Sin embargo, al acoplar este sensor a un microcontrolador que ejecuta un modelo matemático del sistema, podemos expandir sus capacidades.
) solo nos daría información puntual de esa posición y sería un error asumir que toda la barra se encuentra a la misma temperatura. Sin embargo, al acoplar este sensor a un microcontrolador que ejecuta un modelo matemático del sistema, podemos expandir sus capacidades.El objetivo es construir un Sensor Virtual (como el mostrado en la Figura 1) que, utilizando únicamente la lectura ruidosa en y el conocimiento de la física de difusión (matriz F), sea capaz de reconstruir el perfil de temperatura completo a lo largo de la barra, que por practicidad, únicamente se mostraran los nodos marcados.
y el conocimiento de la física de difusión (matriz F), sea capaz de reconstruir el perfil de temperatura completo a lo largo de la barra, que por practicidad, únicamente se mostraran los nodos marcados.
Figura 1: Esquema del problema de estimación de temperatura. Se muestra la evolución libre de una distribución de calor disipándose hacia los extremos fríos. El desafío consiste en estimar la temperatura en todos los puntos deseados (en el problema la barra se divide en 20 puntos con los nodos indicados a graficar) utilizando exclusivamente la lectura del sensor RTD ubicado en la frontera izquierda.
Como se trabajo en la práctica llamada "Filtro de Kalman", este filtro opera en en pasos recursivos en cada instante k [2]. A continuación se introduce nuevamente con una pequeña modificación (tomar en cuenta ya que se deberá responder al finalizar la sección):
% Tuning
Q = 1 * eye(N); % Ruido de proceso
R = sigma_sensor^2; % Ruido de medición
% Inicialización
x_est = zeros(N, steps);
P = 1000 * eye(N); % Alta incertidumbre inicial
primera_medicion = mean(x0);%y_medida(1);
x_est(:, 1) = primera_medicion*ones(N, 1);
% Lazo de Estimación
for k = 2:steps
% 1. Predicción
x_pred = F * x_est(:,k-1);
P_pred = F * P * F' + Q;
% 2. Corrección
S = H * P_pred * H' + R;
K = P_pred * H' / S;
z = y_medida(k); % Dato del sensor en x=0
x_est(:,k) = x_pred + K * (z - H * x_pred);
P = (eye(N) - K * H) * P_pred;
end
indices_interes = [1, 5, 6, 10, 16, 20];
posiciones_fisicas = x_pos(indices_interes);
limite_y_inf = -2;
limite_y_sup = T_inicial + 5;
for k = 1:length(indices_interes)
idx = indices_interes(k);
pos_m = posiciones_fisicas(k);
figure('Name', sprintf('Nodo %d', idx), 'Color', 'w');
plot(t, x_true(idx,:), 'k', 'LineWidth', 2); hold on;
plot(t, x_est(idx,:), 'b--', 'LineWidth', 2);
if idx == 1
plot(t, y_medida, 'r.', 'MarkerSize', 6);
lgd = legend('Temperatura Real', 'Estimación Kalman', 'Medición Ruidosa', ...
'Location', 'Best', 'FontSize', 10);
title_str = sprintf('Nodo 1 (x = %.2f m) - UBICACIÓN DEL SENSOR', pos_m);
else
lgd = legend('Temperatura Real', 'Estimación Kalman', ...
'Location', 'Best', 'FontSize', 10);
title_str = sprintf('Nodo %d (x = %.2f m) - SENSOR VIRTUAL', idx, pos_m);
end
ax = gca;
ax.Color = 'w';
ax.XColor = 'k';
ax.YColor = 'k';
ax.GridColor = 'k';
ax.GridAlpha = 0.15;
title(title_str, 'FontSize', 12, 'FontWeight', 'bold', 'Color', 'k');
ylabel('Temperatura (°C)', 'FontSize', 11, 'Color', 'k');
xlabel('Tiempo (s)', 'FontSize', 11, 'Color', 'k');
set(lgd, 'TextColor', 'k', 'Color', 'w', 'EdgeColor', 'k');
grid on;
ylim([limite_y_inf limite_y_sup]);
xlim([0 T_total]);
end
figure('Name', 'Error Extremo', 'Color', 'w');
error_extremo = abs(x_true(end,:) - x_est(end,:));
area(t, error_extremo, 'FaceColor', [0.85 0.33 0.1], 'FaceAlpha', 0.4, 'EdgeColor', 'r');
ax = gca;
ax.Color = 'w';
ax.XColor = 'k';
ax.YColor = 'k';
ax.GridColor = 'k';
ax.GridAlpha = 0.15;
title('Error Absoluto de Estimación en el Extremo (x=L)', 'FontSize', 12, 'FontWeight', 'bold', 'Color', 'k');
ylabel('|Error| (°C)', 'FontSize', 11, 'Color', 'k');
xlabel('Tiempo (s)', 'FontSize', 11, 'Color', 'k');
grid on;
xlim([0 T_total]);
Actividad
- Cambia el valor de R a un valor muy alto. ¿Qué le pasa a la línea azul (estimación) en el Nodo 1? ¿Se vuelve más suave o más ruidosa? ¿Tarda más en reaccionar?
- Cambia el valor de Q a un valor muy bajo. ¿El filtro sigue al sensor o se queda pegado a lo que dice la física (la difusión suave)?
- Cambia la línea de inicialización para que inicie en Ceros en lugar de usar el promedio o la medición.Mira la gráfica del Nodo 20 (x=1m). ¿Cuánto tiempo tarda el sensor virtual en darse cuenta de que la barra no estaba fría? Esto demuestra la importancia de una buena condición inicial en sistemas lentos.
- Utilice su código original. Modifique la variable L a los siguientes valores: 0.5 m, 2.0 m, y 5.0 m. Para cada caso, ejecute la simulación y observe la gráfica del Nodo 20 (x=L). ¿La estimación (línea azul) converge más rápido o más lento que en el caso original de 1m? ¿Por qué?.
- Regrese la longitud a la incial, y cambie el orden del sistema haciendolo más pequeño y más grande. Al cambiar el orden debe ajustar el código. ¿Cómo cambia la estimación con el orden? ¿Vale la pena aumentar el orden del sistema?
- Cambia el modelo y el filtro de Kalman, de tal manera, que ahora no decaiga a 0, sino a un valor fijado de temperatura, es decir, si la
del pulso corresponde a ese valor, que la temperatura en la barra decaiga a la temperatura ambiente fijada. Nota: la modificación del filtro de Kalman, involucra una entrada.
Reflexión Final
Aunque solo medimos la temperatura en la orilla, el filtro de Kalman fue capaz de reconstruir la temperatura en el centro y el extremo opuesto. Esto es posible porque el modelo matemático le dice al filtro que el calor no puede aparecer mágicamente en la orilla; tuvo que haber viajado desde el centro.
El filtro usa esa correlación física para "invertir" la causa y el efecto: al ver el efecto en la orilla, deduce la causa en el centro, filtrando el ruido en el proceso. Esto convierte a un sensor barato y un microcontrolador en un sistema de instrumentación avanzado.
Referencias
[1] F. P. Incropera y D. P. DeWitt, Fundamentos de transferencia de calor. México: Pearson Educación, 1999
[2] A. D. Núñez-Paredes, R. L. Abad-Rodríguez, R. G. Ramírez-Chavarría and L. Torres, "State-Estimation as a Virtual Sensing Scheme in a Heat Conduction Experiment," en Proceeding of the XX Latin American Control Congress (LACC 2025), Cancún, México, Oct. 15-17, 2025: Springer, en prensa.
[3] R. Haberman, Applied Partial Differential Equations with Fourier Series and Boundary Value Problems, 5ta ed. Upper Saddle River, NJ: Pearson, 2012.
[4] R. J. LeVeque, Finite Difference Methods for Ordinary and Partial Differential Equations: Steady-State and Time-Dependent Problems. Philadelphia, PA: SIAM, 2007.