Práctica 3: Sensores Inteligentes

Estimación de paramétros en línea

• Detectar y compensar deriva o envejecimiento.
• Transformar parámetros estimados en la magnitud de interés (transducción basada en modelo).
• Implementar técnicas de autocalibración y diagnósticos en línea.

Paso 1. Modelo del sensor

Un sensor capacitivo de humedad relativa funciona midiendo los cambios en la capacitancia de un condensador cuya constante dieléctrica varía con el contenido de humedad del aire. El sensor está formado por dos placas conductoras separadas por un material dieléctrico (polímero o cerámica porosa) que absorbe moléculas de agua del ambiente. Cuando la humedad relativa (HR) del aire aumenta, más moléculas de agua se adsorben en el dieléctrico, incrementando su constante dieléctrica . Sin embargo, dado que el capacitor no es perfecto y presenta una componente de disipación, también existe una resistencia eléctrica asociada a él.

La Figura 2 muestra el circuito equivalente del sensor de HR, compuesto por un resistor y un capacitor. El sensor es excitado con un voltaje , y la respuesta del sensor o salida, corresponde al voltaje del capacitor amplificado por un factor de diez (10x).

donde es el voltaje del resistor y es el voltaje del capacitor. Por otro lado, la corriente a través del capacitor está dada por

El coeficiente de la primera derivada en (4) puede definirse como la constante de tiempo del circuito, . Asimismo, debemos recordar que el voltaje del capacitor se amplifica por un factor K, por lo que, la salida del sistema completo es . Con estas consdieraciones, el modelo dinámico puede escribirse finalmente en terminos del voltaje de salida como

• A partir de la ecuación diferencial (5) obtenga su función de transferencia.

Paso 2. Simulación del sensor

clear; clc; close all;

 

%% Parámetros del circuito

R = 10e3; % Resistencia [ohms]

C = 100e-9; % Capacitancia [farads]

tau = R * C; % Constante de tiempo [s]

K = 10; %ganancia de la etapa de amplificación 10x

fprintf('Constante de tiempo τ = %.3e s\n', tau);

Constante de tiempo τ = 1.000e-03 s

 

%% Funciones de transferencia

s = tf('s');

 

% Salida = voltaje en el capacitor

G_C = K / (tau*s + 1);

 

 

%% Respuesta temporal a un escalón de 5 V

Ts = 0.0001; %Periodo de muestreo

t = 0:Ts:0.025; % tiempo de simulación [s]

N = length(t);

Vi = 0.5*sin(2*pi*100*t);%0.01*ones(size(t)); % entrada: escalón de 5 V

 

% Simulación del sistema

A = -1/tau;

B = K/tau;

Vo = zeros(size(t));

for k = 2:length(t)

Vo(k) = (Vo(k-1) + Ts*(A*Vo(k-1) + B*Vi(k-1)));

end

 

r = 1e-1*randn(1,N); %ruido de medición

 

Vo = Vo + r; % Voltaje de salida medido

 

 

figure('Color', 'w');

plot(t, Vi, 'b--', 'LineWidth', 2);

hold on;

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

 

set(gca, 'Color', 'w', ...

'XColor', 'k', ...

'YColor', 'k', ...

'GridColor', 'k', ...

'GridAlpha', 0.2);

 

xlabel('Tiempo [s]', 'Color', 'k');

ylabel('Voltaje [V]', 'Color', 'k');

title('Respuesta temporal del circuito RC', 'Color', 'k');

lgd = legend('Entrada V_{i}(t)', 'Voltaje en el capacitor V_o(t)', 'Location', 'southeast');

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

grid on;

hold off;

 

• Modifica el código anterior para que ahora el voltaje de entrada sea una señal senoidal de amplitud 0.5 y frecuencia de 100 Hz. El ruido de medición deberá tener una amplitud de 0.1.
• ¿Cómo es el voltaje de salida respecto al voltaje de entrada?

Paso 3. Problema de estimación

donde y son los parámetros a estimar. La función de transferencia en (6) puede reescribirse en términos de la señal de salida como

lo cual corresponde a la transformada de Laplace de la ODE en (5). Comparando (5) y (7) es claro observar que el término corresponde a la derivada del voltaje de salida. Para obtener dicha señal, es necesario realizar la diferenciación numérica de ; sin embargo, en caso de emplear métodos numéricos convencionales, ésto conlleva a amplificar el ruido de medición, lo que resulta poco viable en la práctica.

Para obtener la derivada del voltaje, podemos emplear un filtro lineal de primer orden con frecuencia de corte , cuya función de transferencia es

donde , y , son las señales filtradas del voltaje de entrada, de salida y de la derivada del voltaje de salida, respectivamente. De forma interesante, el modelo en la ecuación (9) es lineal en los parámetros b y a. De esta forma, (9) puede escribirse en notación matricial como

donde se denomina regresor, y es el vector de parámetros a estimar (desconocido), con ⊤ la transpuesta del vector. Con ello, el problema de estimación de parámetros del circuito equivalente del sensor de HR, puede escribirse como un problema de regresión lineal de la forma

%% Filtrando las señales de entrada y salida

w_c = 2.5; % polo del filtro

 

 

Vof = zeros(size(t));

Vif = zeros(size(t));

Vof_dot = zeros(size(t));

 

% Se aplica el filtro muestra por muestra de ambas señales

% Esto es equivalente a filtrar en tiempo real

for k = 2:length(t)

Vof(k) = Vof(k-1) + Ts*(-w_c*Vof(k-1) + w_c*Vo(k));

Vif(k) = Vif(k-1) + Ts*(-w_c*Vif(k-1) + w_c*Vi(k));

 

% Aproximación de la derivada dVo/dt (medida) usando diferencias

% finitas

Vof_dot(k) = (Vof(k) - Vof(k-1))/Ts;

 

end

 

 

 

figure('Color', 'w');

subplot(3,1,1)

plot(t, Vof_dot, '-k', 'LineWidth', 1.2)

grid on

set(gca, 'Color', 'w', 'XColor', 'k', 'YColor', 'k', 'GridColor', 'k');

xlabel('Tiempo [s]', 'Color', 'k');

ylabel('Voltaje [V]', 'Color', 'k');

title('Derivada de la señal de salida filtrada', 'Color', 'k');

subplot(3,1,2)

plot(t, Vof, '--b', 'LineWidth', 1.2)

grid on

set(gca, 'Color', 'w', 'XColor', 'k', 'YColor', 'k', 'GridColor', 'k');

xlabel('Tiempo [s]', 'Color', 'k');

ylabel('Voltaje [V]', 'Color', 'k');

title('Señal de salida filtrada', 'Color', 'k');

subplot(3,1,3)

plot(t, Vif, '--r', 'LineWidth', 1.2)

grid on

set(gca, 'Color', 'w', 'XColor', 'k', 'YColor', 'k', 'GridColor', 'k');

xlabel('Tiempo [s]', 'Color', 'k');

ylabel('Voltaje [V]', 'Color', 'k');

title('Señal de entrada filtrada', 'Color', 'k');

 

% Construcción del regresor

Phi = [Vif' -Vof'];

 

%Solución del problema de regresion por mínimos cuadrados ordinarios

theta_LS = pinv(Phi'*Phi)*Phi'*Vof_dot';

fprintf('Parámetro b estimado = %.2f \n', theta_LS(1));

Parámetro b estimado = 9915.60

fprintf('Parámetro a estimado = %.2f \n', theta_LS(2));

Parámetro a estimado = 994.03

• Investiga y deduce la ecuación normal del estimador de mínimos cuadrados para el problema de regresión
• A partir de los parametros estimados en el código anterior, calcula el error entre estos y los parámetros verdaderos. Recuerda que los parámetros verdaderos están dados por y .

Paso 3. Mínimos cuadrado recursivos

El algorimto trabaja con la ésima muestra de las señales , y , una matriz de covarianza inicial , el vector de parámetros estimados inicial , y un hiperparámetro llamado factor de olvido. De esta forma el RLS está conformado por los siguientes seis pasos:

1. Construir el regresor:
2. Cálculo del término de corrección:
3. Cálculo de la señal de salida estimada:
4. Cálculo del error de estimación:
5. Actualización de los parámetros:
6. Actualización de la matriz de covarianza:

los cuales se realizan de forma iterative desde hasta N, donde Nes el número total de muestras de las señales , y .

% Algoritmo RLS

 

% Parámetros "verdaderos"

a_true = 1/tau;

b_true = K/tau;

 

 

% Inicialización

theta = [0; 0]; % Estimación inicial [b_hat; a_hat]

P = 1E6*eye(2); % Covarianza inicial

lambda = 0.95; % Factor de olvido

 

theta_hist = zeros(2, length(t));

 

 

for k = 2:N

 

phi = [Vif(k); -Vof(k)];

Kappa = P*phi/(lambda + phi'*P*phi);

Vo_hat = phi'*theta;

err = Vof_dot(k) - Vo_hat;

theta = theta + Kappa*err;

P = (P - Kappa*phi'*P)/lambda;

theta_hist(:,k) = theta; %almacena la evolución de los parámetros estimados

end

 

%% Grafica resultados

figure('Color', 'w');

subplot(2,1,1);

plot(t, theta_hist(1,:), 'r', 'LineWidth', 1.5); hold on;

yline(b_true, '--b', 'LineWidth', 1.5);

set(gca, 'Color', 'w', 'XColor', 'k', 'YColor', 'k', 'GridColor', 'k');

title('Parámetro b estimado', 'Color', 'k');

lgd1 = legend('b_{est}', 'b_{true}');

set(lgd1, 'TextColor', 'k', 'Color', 'w', 'EdgeColor', 'k');

ylim([min(theta_hist(1,:))-1000, max(theta_hist(1,:)) + 1000])

 

grid on;

 

subplot(2,1,2);

plot(t, theta_hist(2,:), 'r', 'LineWidth', 1.5); hold on;

yline(a_true, '--b', 'LineWidth', 1.5);

set(gca, 'Color', 'w', 'XColor', 'k', 'YColor', 'k', 'GridColor', 'k');

title('Parámetro a estimado', 'Color', 'k');

lgd2 = legend('a_{est}', 'a_{true}');

set(lgd2, 'TextColor', 'k', 'Color', 'w', 'EdgeColor', 'k');

ylim([min(theta_hist(2,:))-1000, max(theta_hist(2,:)) + 1000])

grid on;

El resultado del código anterior es la señal de evolución en el tiempo de los parámetros estimados del modelo y .

• Identifica los pasos del algoritmo RLS en el código anterior y escribe los comentarios pertinentes para describir su funcionamiento.
• ¿Cómo es el valor de los parámetros estimados respecto al de los parámetros reales? ¿Hay alguna forma de cuantificar esa diferencia?
• Modifica el código anterior para cambiar el valor de lambda a 0.1, 0.5 y 1.0. ¿Cuál es el efecto de lambda en la estimación?
• Modifica el código anterior para cambiar el valor de la matriz de convariance a 1, 1e3 y 1e5. ¿Cuál es el efecto de P en la estimación?

Paso 4: Interpretación de los parámetros

Lo anterior, corresponde a un sistema de dos ecuaciones con dos incógnitas, las cuales pueden ser obtenidas a partir de los parámetros estimados y . De esta forma se tiene que la constante de tiempo del sensor puede estimarse como

• Modifica el código de la sección anterior considerando P = 1E6*eye(2) y lambda = 0.95 calcular los valores de y como función del tiempo. Es decir, estos deberá evolucionar en el tiempo como una señal, similar a lo obtenido en los parámetros .

Reflexión final

Ejercicio para el estudiante

donde el superíndice ⋆ hace referencia a que b y a no son constantes en todo el intervalo de tiempo. Físicamente, esto implicará que desde un inicio y hasta los 0.0125 segundos, los parámetros b y a dependerán de los valores ms y 10. A partir de 0.0126 segundos, la capacitancia cambió debido a un cambio en la HR, y la constante de tiempo debe ser ms, mientras que la ganancia también cambio a un tercio de su valor origina. El siguiente código te ayudara a definir dichos cambios en los parámetros del modelo.

% ****************************************************************

% 1. Simulación del sistema

% ****************************************************************

 

Vo = zeros(size(t));

for k = 2:length(t)

 

%Vo(k) = (Vo(k-1) + Ts*(A*Vo(k-1) + B*Vi(k-1)));

if k >round(N/2)

A = -1/(2*tau);

B = (K/3)/tau;

end

 

Vo(k) = (Vo(k-1) + Ts*(A*Vo(k-1) + B*Vi(k-1)));

end

 

r = 1e-1*randn(1,N); %ruido de medición

 

Vo = Vo + r; % Voltaje de salida medido

 

figure('Color', 'w');

plot(t, Vi, 'b--', 'LineWidth', 1.2);

hold on;

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

set(gca, 'Color', 'w', ...

'XColor', 'k', ...

'YColor', 'k', ...

'GridColor', 'k', ...

'GridAlpha', 0.2);

xlabel('Tiempo [s]', 'Color', 'k');

ylabel('Voltaje [V]', 'Color', 'k');

title('Respuesta temporal del circuito RC', 'Color', 'k');

lgd = legend('Entrada V_{i}(t)', 'Voltaje en el capacitor V_o(t)', 'Location', 'southeast');

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

 

grid on;

hold off;

 

% ****************************************************************

% 2. Se aplica el filtro muestra por muestra de ambas señales

 

% Aquí va el código para filtrar las señales

% ****************************************************************

 

 

% ****************************************************************

% Algoritmo RLS

 

% Aquí va el código del estimador RLS

 

% ****************************************************************

Actividad

1. Compara los resultados del algoritmo RLS para lambda 0.1, 0.5 y 0.995. ¿Cuál es la diferencia entre los resultados obenitos?
2. Cambia la señal de entrada por un escalón de amplitud 0.5 volts. ¿Qué sucedió con el desemepeño del algoritmo RLS?
3. Investigue algunas aplicaciones del algoritmo RLS en áreas como procesamiento digital de señales y control automático.

Referencias