Práctica 2: Sensores Inteligentes

Sensor Inteligente de Desplazamiento

*Objetivo:* Comprender el comportamiento dinámico de un sensor mecánico modelado como un sistema masa–resorte–amortiguador (MRA), identificando la influencia de los parámetros físicos (masa, constante de resorte y coeficiente de amortiguamiento) sobre la respuesta temporal y frecuencial del sistema, y su utilidad como mecánismo de transducción en sensores inteligentes.

*Modelo dinámico del sensor:* De acuerdo con los primero principios que rígen la física de un sistema MRA, el modelo de éste está dado por la siguiente ecuación diferencial ordinaria (ODE: ordinary differential equation) [1]:

La ODE anterior permite describir como evolucionan las variables de posición

, velocidad

, y aceleración

, de acuerdo con los parámetros masa m, coeficiente amortiguamiento b y la constante del resorte k, así como con la fuerza aplicada

*Motivación:* En esta práctica estudiaremos el efecto de los parámetros del modelo dinámico en la respuesta del sensor. Con ello, será posible establecer un marco de referencia para comprender el funcionamiento de sensores inteligentes basados en modelo [2].

*Procedimiento:* La práctica está dividida en seis pasos. Lee detenidamente la descripción de cada paso, y en caso de que exista un fragmento de código, deberás transcribirlo en un archivo .mlx (LiveScript) de MATLAB, y responder las preguntas de cada sección.

Paso 1. Generación de datos

Como primer paso, emplearemos las medidas de desplazamiento o posición obtenidas con un sensor cuya respuesta se aproxima a la de un sistema MRA. La configuración experimental para obtener las señales de desplazamiento se muestra en la Figura 1.

El siguiente fragmento de código se encarga de cargar los datos del experimento, almacenados en una matriz cuya primera columna corresponde al tiempo y la segunda a los datos de la posición. Copia y pega el siguiente fragmento de código en un LiveScript de MATLAB con extensión .mlx

clear all;

close all;

clc

datos = load("masa_resorte_tracker.mat"); %lectura de datos

t = datos.t_aire; % Columna de datos de tiempo

Ts = t(2) - t(1); % Periodo de muestreo

y= datos.y_aire; % Columna de datos de posición

y = detrend(y); % Pre-procesamiento de la señal p/quitar nivel DC

figure('Color', 'w');

plot(t, y, 'r', 'LineWidth', 1.5);

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

'XColor', 'k', ...

'YColor', 'k', ...

'GridColor', 'k');

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

ylabel('Desplazamiento [m]', 'Color', 'k')

title('Señal medida por el sensor', 'Color', 'k')

grid on

Preguntas de la sección:

¿Qué patrón observa en la señal de desplazamiento (posición)?
¿Cuál es la duración de dicha señal en segundos?
¿Cuál es el valor RMS del desplazamiento medido?

Paso 2. Estructura del modelo

Posteriormente, es necesario establecer una estructura del modelo MRA en variables de estado [1]. Para ello, se definen las variables de estado como

Por simplicidad usaremos la notación

, para la primera derivada y

para la segunda derivada. Así, definimos el vector de estado como

La idea es entonces a la ODE del sistema MRA en (1) como un modelo en espacio de estados (SS: state-space) para un sistema lineal e invariante con el tiempo (LTI: linear time-invariant)

Sustituyendo la definición de las variables de estado en (2) en la ecuación diferencial (1), y arreglando los términos en notación matricial, la ecuación de estado del modelo es

Mientras que la ecuación de salida corresponde a la medición de la posición, y está dada por

Preguntas de la sección:

¿Cómo se modificaría la ecuación (6) sí la variable a medir fuera la velocidad?
¿Cómo se obtendría la aceleración a partir de la velocidad?

Paso 3. Definir el problema de estimación

3.1. Estimación del modelo

Una vez establecida la estructura del modelo dinámico, es necesario plantear el problema de estimación del modelo. Esto es, buscaremos los parámetros de la ecuación de estado (5), a partir de la señal de posición medida experimentalmente. Para ello, haremos uso de un método de identificación basado en subespacios denominado N4SID (numerical algorithm for state space subspace system identification).

Por simplicidad, no profundizaremos en el fundamento teórico de N4SID, pero se deja al lector la posibilidad de consultar la bibliografía [3] para mayor referencia.

En el ambiente MATLAB, el System Identification Toolbox [4], es capaz de realizar la identificación del modelo dado por (5) y (6) a través del siguiente fragmento de código, el cual deberás trasncribir en tu archivo .mlx

nx = 2;    %Orden del modelo MRA
data_id = iddata(y,[],Ts); %creación del objeto de identificación a partir de la señal medida 
 
% Opciones para el algoritmo N4SID)
Options = n4sidOptions;                                      
opt.N4weight = 'CVA'; 
opt.N4Horizon = 'auto';
opt.EstimateCovariance = 'true';
 
% Modelo identificado por N4SID
Modelo = n4sid(data_id,nx,'Form','canonical','Ts',0,Options)
Modelo =
  Continuous-time identified state-space model:
      dx/dt = A x(t) + K e(t)
       y(t) = C x(t) + e(t)
 
  A = 
             x1        x2
   x1         0         1
   x2    -13.98  -0.06109
 
  C = 
       x1  x2
   y1   1   0
 
  K = 
          y1
   x1  33.13
   x2  90.61
 
Parameterization:
   CANONICAL form with indices: 2.
   Disturbance component: estimate
   Number of free coefficients: 4
   Use "idssdata", "getpvec", "getcov" for parameters and their uncertainties.

Status:                                             
Estimated using N4SID on time domain data "data_id".
Fit to estimation data: 93.3% (prediction focus)    
FPE: 5.496e-06, MSE: 5.443e-06                      
 
Model Properties

Como resultado, se puede observar que N4SID logró identificar un modelo SS, cuya estructura coincide con la presentada en las ecuaciones (5) y (6).

Preguntas de la sección:

¿De qué orden fue el modelo estimado por N4SID? Escribe en tu archivo .mlx en formato de ecuación, las matrices A y C estimadas por N4SID.
¿Cuál fue el porcentaje de ajuste de los datos al modelo? Lo consideras, ¿bueno o malo?

3.2. Evaluación del modelo estimado

Para comparar el resultado anterior, graficaremos la señal de posición medida experimentalmente

, con la señal estimada por el modelo

. Para ello, transcribe el siguiente fragmento de código a tu archivo .mlx

kstep=1;

[y_hat,fit_y,ic]=compare(y,Modelo,kstep);

figure('Color', 'w');

subplot(2,1,1)

plot(y, '-', 'LineWidth', 3)

hold on

plot(y_hat, '-', 'LineWidth', 0.7)

hold off

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

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

ylabel('Desplazamiento [m]', 'Color', 'k')

title('Señal de desplazamiento', 'Color', 'k')

lgd = legend('$\mathbf{y}$', '$\mathbf{\hat{y}}$', 'Interpreter', 'Latex');

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

grid on

subplot(2,1,2)

plot(y - y_hat, '-r', 'LineWidth', 2)

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

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

ylabel('Desplazamiento [m]', 'Color', 'k')

title('Error de estimación', 'Color', 'k')

grid on

Preguntas de la sección:

¿Cuál es la función del comando compare()?
Investigue el concepto de error cuadrático medio (RMSE) y calculelo para la señal de error de estimación mostrado en la gráfica.

Paso 4. Análisis de los parámetros estimados

Una vez identificado el modelo SS a partir de los datos experimentales, realizaremos el análisis de los parámetros estimados en la matriz del sistema. De acuerdo con el resultado, el objeto Modelo, muestra que la matriz estimada del sistema tiene la siguiente estructura

Comparando la expresión (7) con la ecuación del sistema en (5), se obtiene directamente que

Lo que resulta en un sistema de dos ecuaciones con tres incógnitas. Sin embargo, previo al experimento el valor de la masa fue medido con un valor de 0.6984 kg. Al sutituir dicho valor en (8) y (9), es entonces posible calcular los valores de la constante del resorte, k y del coeficiente de amortiguamiento b. Para ello, es importante que escribas nuevamente la matriz

estimada en el paso 3.

El siguiente fragmento de código nos ayudara a extraer los valores numéricos de

[A_hat,~,C_hat,~,K] = idssdata(Modelo);  % Extrae la matriz del sistema A
a1_hat = A_hat(2,1);                     % Extrae el elemento (2,1) de la matriz asociado con el parámetro -k/m
a2_hat = A_hat(2,2);                     % Extrae el elemento (2,2) de la matriz asociado con el parámetro -b/m

Preguntas de la sección:

Con los valores numéricos obtenidos del código anterior, y sabiendo que la masa es de 0.6984 kg. Escriba el código necesario para obtener los valores de k y m a partir de las ecuaciones (8) y (9). Utilice el comando fprintf() para mostrar sus valores en la consola.

Paso 5. Analizar el efecto del amortiguamiento

Para observar el efecto del amortiguamiento estimado b, podemos realizar un análisis del modelo SS identificado en el dominio de la frecuencia. Para ello, es posible emplear la representación en función de transferencia del sistema MRA, la cual tiene la siguiente estructura

Usando los valores de la masa conocida y los valores estimados de la constante del resorte y el coeficiente de amortiguamiento obtenidos en el paso anterior, en la función de transferencia (10), podemos realizar el análisis en frecuencia con el siguiente fragmento de código.

s = tf('s'); m = 0.6984;
G = 1/(m*s^2 + b_hat*s + k_hat);     % Función de transferencia
Unrecognized function or variable 'b_hat'.
% Al ejecutar el código, un error aparecerá
% deberá obtener los parámetros físicos señalados e introducirlos en el programa.
 
figure
bode(G)                              % Respuesta en frecuencia
grid on
title('Respuesta en frecuencia del sistema masa-resorte-amortiguador')
 
wn = sqrt(k_hat/m);                 % Estimación de la frecuencia natural
fprintf('Frecuencia natural ω_n = %.2f rad/s\n', wn);

Como resultado podemos obtener un parámetro más, llamado frecuencia natural

, el cuál servirá como una característica adicional que el sensor puede extraer.

Paso 6. Reflexión final

En esta práctica estudiamos como procesar, analizar y trabajar la señal de un sensor de desplazamiento para otorgarle inteligencia en el sentido de extraer la mayor cantidad de información posible a partir de una sola señal medida experimentalmente. Con ella, fue posible identificar un modelo dinámico de un sistema MRA, cuyos parámetros estimados nos dieron más información que la señal por sí sola. Posteriormente, con el modelo identificado también pudimos encontrar un parámetro adicional en el dominio de la frecuencia. Así, la señal del sensor pudo proveer, mediante algoritmos de estimación, un conjunto de tres parámetros

los cuales son características únicas de la medición que el sensor está realizando, y pueden servir como mecanismo de clasificación de muestras, por ejemplo.

Ejercicio para el estudiante

Repita el procedimiento anterior, de los puntos 1 al 6, pero ahora considerando que la masa del sistema se desplaza en agua. Observe que la diferencia entre el experimento anterior y este, es el medio en el que la masa se desplaza. La hipótesis a verificar es como las propiedades físicas del medio determinarán los parámetros del modelo MRA, y cómo éstos pueden servir como mecánismo de transducción.

El código para realizar la lectura de datos de desplazamiento en agua es el siguiente.

datos = load("masa_resorte_tracker.mat"); %lectura de datos 
 
t = datos.t_agua;     % Columna de datos de tiempo
Ts =  t(2) - t(1);    % Periodo de muestreo
y= datos.y_agua;      % Columna de datos de posición en aire
y = detrend(y);       % Pre-procesamiento de la señal p/quitar nivel DC

Actividad

Estime el modelo en espacio de estados usando el algorimto N4SID.
Obtenga los parámetros para las mediciones en aire, y comparelos con los obtenidos en agua.
Obtenga la respuesta en frecuencia del modelo en aire, y compare la frecuencia de amortiguamiento estimada.
¿Cómo afecta el coeficiente de amortiguamiento al tiempo de estabilización del sistema?
¿Qué parámetro del modelo modifica la frecuencia natural del sistema?
¿Cómo se relaciona el parámetro anterior con las propiedades físicas del medio?

Referencias

[1] Ogata, Katsuhiko. Ingeniería de control moderna (5.ª edición). Pearson Educación, 2010.

[2] Bhuyan, M. (2011). Intelligent Instrumentation: Principles and Applications. Boca Raton, FL: CRC Press. ISBN 978-1-4200-8953-0.

[3] Principles of System Identification: Theory and Practice — Arun K. Tangirala, CRC Press, 2015.

[4] System Identification Toolbox for Use with MATLAB: User's Guide — Lennart Ljung, MathWorks, 2004.