Notas de clase: Modelación experimental

Table of Contents

Notas de clase: Modelación experimental Clase 1: Presentación del curso Tema: Objeto y método de la modelación experimental Clase 2: Objeto y método de la modelación experimental Clase 3: Objeto y método de la modelación experimental Clase 4: Objeto y método de la modelación experimental Estructuras de modelos Tema: Métodos no paramétricos de identificación Clase 5: Métodos no paramétricos de la respuesta temporal Respuesta temporal para sistemas de primer orden Ajuste mediante modelos de primer orden Respuesta temporal para sistemas de segundo orden Ajuste mediante modelos de segundo orden Resumen de los métodos de respuesta temporal Clase 6: Método no paramétrico de la correlación Estimación de la secuencia de ponderación mediante correlación Tema: Métodos paramétricos de identificación Clase 7: Método de mínimos cuadrados (teoría) Clase 8: Método de mínimos cuadrados (aplicaciones) Clase 9: Métodos de variable instrumental Clase 10: Métodos de variable instrumental Funciones adicionales para el curso

Clase 1: Presentación del curso

Mapa mental con las ideas de la clase.

Tema: Objeto y método de la modelación experimental

Clase 2: Objeto y método de la modelación experimental

Clase 3: Objeto y método de la modelación experimental

Clase 4: Objeto y método de la modelación experimental

Estructuras de modelos

OE: No se desea modelar la perturbación, se considera que ésta es un ruido que se adiciona a la salida. nk es el retardo puro del sistema.

ARX (auto-regresivo con variable exógena): Es la estructura más común (se relaciona con el método de mínimos cuadrados), útil si se tiene una buena relación ruido/señal. La perturbación no se modela pero actúa por el canal de entrada del sistema. nk es el retardo puro del sistema.

ARMAX (auto-regresivo, media móvil con variable exógena: útil si la perturbación actúa por el canal de entrada (perturbaciones en la carga).

Box-Jenkins:

Estructura general:

Tema: Métodos no paramétricos de identificación

Clase 5: Métodos no paramétricos de la respuesta temporal

Respuesta temporal para sistemas de primer orden

Una función de transferencia para un sistema de primer orden tiene la forma

, siendo

. Dicha función de transferencia corresponde entonces con la ecuación diferencial

cuya solución viene dada por

para el caso de una entrada escalón con amplitud A,

En dicho contexto, es posible dar a los parámetros de la función de transferencia una interpretación acorde con la forma de la solución de la ecuación diferencial:

constante de tiempo del sistema: se deduce haciendo , en cuyo caso . De forma tal, T es el tiempo que tarda el sistema en ir de su valor inicial al de su valor final (de estabilización).
ganancia del sistema: se deduce calculando el . De forma tal, k determina la amplificación con la que el sistema responde a la entrada.
retardo puro del sistema
amplitud de la entrada escalón aplicada al sistema

Ejemplo: comparación de sistemas de primer orden

¿Qué podría decirse de los siguientes sistemas de primer orden?

Comparando en MatLab: con la función tf podemos crear funciones de transferencia especificando los coeficientes del polinimio del numerador y del denominador.

clear;clc;close all
tf1=tf([0.5],[1 0.25])
tf1 =
 
    0.5
  --------
  s + 0.25
 
Continuous-time transfer function.
tf2=tf(2,[1 4])
tf2 =
 
    2
  -----
  s + 4
 
Continuous-time transfer function.

Con la función step podemos graficar rápidamente la respuesta del sistema a una entrada de tipo escalón unitario.

figure(1)

clf

step(tf1,20)

hold on

step(tf2,20)

grid on

xline([4,0.25],'-',{'T_1=4','T_2=0.25'})

legend({'G_1(s)','G_2(s)'})

En conclusión, el primer sistema crece más lento pero se estabiliza en un mayor valor mientras que el segundo crece más rápido y se estabiliza en un menor valor.

Ajuste mediante modelos de primer orden

Método visual: Los valores de los 3 parámetros se deducen de la gráfica

Método de regresión lineal: Los valores de los 3 parámetros se deducen mediante regresión. Se parte de la estructura

y se deduce que

, o visto de otra forma, el problema se ha llevado a la forma

. El valor de k se calcula fácilmente como

Ejemplo:

Se genera un modelo de segundo orden especificando los ceros, polos, ganancia y retardo del sistema mediante la función zpk. El sistema que elegimos para el ejemplo no tiene ceros, tiene dos polos (-1.8 y -2.2), tiene una ganancia de 1.24 y un retardo de 0.89 segundos.

%creación de un modelo de segundo orden

G=zpk([],[-1.8 -2.2],1.24,'InputDelay',0.89)

G =
 
                      1.24
  exp(-0.89*s) * ---------------
                 (s+1.8) (s+2.2)
 
Continuous-time zero/pole/gain model.

%generación de la respuesta temporal

figure(1)

clf

step(G);%comprobamos que el sistema se estabiliza

%[y,t]=step(G);%recolectamos los datos tal y como los arroja matlab

[y,t]=step(G,[0:0.01:6]);%recolectamos los datos a nuestro gusto

inicio=find(y>0,1)%identificamos el dato a partir del cual el modelo responde a la entrada

inicio = 91

yvals=y(inicio:end);%valores útiles de y

tvals=t(inicio:end);%valores útiles de t

k=yvals(end)-y(1);%valor de k según la fórmula

Método de regresión

[p,s]=polyfit(tvals,log(1-yvals/(k+eps)),1);%regresión lineal al logaritmo de y
T=-1/p(1);tau=-p(2)/p(1);%valores de parámetros  T y tau según las fórmulas
tf_est=tf(k,[T 1],'InputDelay',tau)%función de transferencia resultante
tf_est =
 
                    0.313
  exp(-1.44*s) * ------------
                 0.5551 s + 1
 
Continuous-time transfer function.

Método gráfico

tau2=t(inicio);%primer valor donde y reacciona a la entrada
valT=find(y>0.632*y(end),1);%posición donde y alcanza el 63.2% de su valor final
T2=t(valT)-tau2;%valor de T por definición
tf_est2=tf(k,[T2 1],'InputDelay',tau2)%función de transferencia resultante
tf_est2 =
 
                  0.313
  exp(-0.9*s) * ----------
                1.08 s + 1
 
Continuous-time transfer function.

Comparación entre métodos

figure(1)

hold on

step(tf_est)

step(tf_est2)

legend({'Original','Regresión','Gráfico'},'Location','best')

Respuesta temporal para sistemas de segundo orden

Los sistemas de segundo orden son de la forma

. Dicha función de transferencia se corresponde con la ecuación diferencial

El caso que nos interesa en esta ecuación ocurre cuando

. ya que las raíces del polinomio característico serán imaginarias y la solución de la ecuación diferencial vendría dada por

De forma similar al caso de la respuesta temporal en sistemas de primer orden, es posible dar interpretación a los parámetros del modelo de orden 2 como sigue

ganancia del sistema
retardo puro del sistema
frecuencia natural (no amortiguada) del sistema
frecuencia amortiguada del sistema
razón de amortiguamiento
amplitud de la entrada escalón del sistema

Dependiendo de los valores de ζ tendremos 4 casos:

, sistema no amortiguado
, sistema subamortiguado
, sistema críticamente amortiguado
, sistema sobreamortiguado

a=4.3;

b=3.1;

c=5.1;

tau=1;

tf1=tf(a,[1 b c],'InputDelay',tau);%función de transferencia a partir de los datos

figure(1)

clf

step(tf1);

Warning: MATLAB has disabled some advanced graphics rendering features by switching to software OpenGL. For more information, click here.

title(["\omega_0= "+sqrt(c)+" k= "+a/c+" \zeta= "+b/(2*sqrt(c))])%cálculo de valores de los parámetros

Ajuste mediante modelos de segundo orden

Para ajustar el modelo de segundo orden utilizamos un método gráfico. Para tal fin, determinamos 4 valores: el retardo del modelo τ, el valor final del modelo

, el sobreimpulso máximo

y el tiempo de pico

[y,t]=step(tf1);%recolección de datos por recomendación de MatLab

inicio=find(y>0,1);%dato a partir del cual el sistema reacciona

tau=t(inicio)%tiempo a partir del cual el sistema reacciona

tau = 1.0102

deltaY=y(end)-y(1)%delta de y por definición

deltaY = 0.8412

Mp=max(y)-y(end);%sobreimpulso máximo pode definición

Tp=t(y==max(y))-tau;%tiempo de pico, se le debe restar el retardo

%zeta=1/sqrt(1+(pi/log(Mp/deltaY))^2)%valor de z (diapositivas)

zeta=-log(Mp)/sqrt(pi^2+log(Mp)^2)%valor de z (Soderstrom)

zeta = 0.7015

omega0=pi/(Tp*sqrt(1-zeta^2))%valor de omega0 por fórmulas

omega0 = 2.3184

tf_est=tf(deltaY*omega0^2,[1 2*zeta*omega0 omega0^2],'InputDelay',tau)%función de transferencia estimada

tf_est =
 
                         4.522
  exp(-1.01*s) * ---------------------
                 s^2 + 3.253 s + 5.375
 
Continuous-time transfer function.

figure(2)

clf

step(tf1)

hold on

step(tf_est)

xline(Tp+tau)

legend({'Original','Estimado'},'Location','best')

Resumen de los métodos de respuesta temporal

Se identifica si el sistema presenta oscilaciones o no. A los modelos que oscilan los aproximaremos por un sistema de orden 2, todos los demás modelos serán aproximados por un modelo de orden 1.
Seleccionar el método para identificar parámetros: gráfico o de regresión.
Implementar las fórmulas para identificar los valores de parámetros.
Validar los resultados de la función de transferencia estimada versus la curva original.

Clase 6: Método no paramétrico de la correlación

Estimación de la secuencia de ponderación mediante correlación

clear;clc;close all

N=1950;%número de muestras

M=10;%Longitud secuencia de ponderación

noise=0.1;%amplitud del ruido

num=[1];%numerador de la tf

den=[1 0.5];%denominador de la tf

tf1=tf(num,den,1)%función de transferencia objetivo

tf1 =
 
     1
  -------
  z + 0.5
 
Sample time: 1 seconds
Discrete-time transfer function.

u=2*(prbs(5,N,[1 0 1 0 0])-0.5);%señal prbs de N datos con valores entre -1 y 1

xdata=0:1:N-1;%instantes de muestreo

[y,t]=lsim(tf1,u,xdata);%simulación del modelo con la entrada diseñada

y=y+noise*rand(length(y),1);%adición del ruido a la entrada

figure(1)

clf

subplot(1,2,1), plot(t,u)%gráfica de la entrada

subplot(1,2,2), plot(t,y)%fráfica de la salida

[cyu,lagsyu] = xcorr(y,u,M-1);%M elementos de la correlación entre y y u

[cu,lagsu] = xcorr(u,0);%correlación de la entrada con ella misma en el instante 0

tf_est=cyu(lagsyu>=0)./(N:-1:N-M+1)'/(cu/N);%fórmula para la secuencia de ponderación

delay=length(den)-length(num);%cálculo del retardo del sistema

tf_real=ldiv([1],[1 0.5],M-delay);%cálculo de la secuencia de ponderación

tf_real=[zeros(1,delay) tf_real];%adición del retardo a la secuencia de ponderación

figure(2)

clf

stem(lagsyu,cyu)%gráfica de la correlación de y y u

title('R_{yu}(\tau)')

pdata=0:1:M-1;

figure(3)

clf

subplot(2,1,1)

stem(pdata,tf_est)%secuencia estimada

title('Estimado')

subplot(2,1,2)

stem(pdata,tf_real)%secuencia real

title('Real')

data = iddata(y,u');%datos para usar las funciones mejoradas de MatLab

sys = impulseest(data,M);%secuencia de ponderación estimada por MatLab

figure(4)

clf

subplot(2,1,1)

h = impulseplot(sys);%secuencia de ponderación estimada por MatLab

subplot(2,1,2)

stem(pdata,tf_real)%secuencia real

title('Real')

figure(5)

clf

response=cumsum(tf_est.*ones(length(tf_est),1));%convolución a la entrada de tipo escalón para la secuencia estimada

stem(pdata,response)

hold on

response2=cumsum(tf_real'.*ones(length(tf_est),1));%convolución a la entrada de tipo escalón para la secuencia real

stem(pdata,response2)

legend({'Secuencia estimada', 'Secuencia real'})

Tema: Métodos paramétricos de identificación

Clase 7: Método de mínimos cuadrados (teoría)

Partimos de un modelo cuya función de transferencia está dada por

Ahora bien, consideraremos que la señal de salida de nuestro sistema estará dada por la superposición de la salida del proceso

y una perturbación estocástica

, es decir

La tarea que nos proponemos entonces es determinar los parámetros desconocidos

tras la recolección de N datos de la entrada

y la salida

del sistema, sujeta a las siguientes suposiciones:

El proceso se encuentra en estado estacionario para .
El orden na y retardo nk del modelo son conocidos exáctamente (la determinación del orden y retardo apropiados se verá más adelante).
La entrada es conocida con precisión.
La perturbación es estacionaria y .

Podemos transformar (1) al dominio del tiempo para obtener la ecuación en diferencias

Adicionalmente, podemos interpretar que (3) es una prediccón de

basada en las mediciones obtenidas hasta el instante anterior, es decir

donde

Pero, por la definición de cada

tenemos que

Siendo

De esta forma, la diferencia entre las observaciones y las predicciones del modelo vendrían dadas por

Donde es de particular importancia el caso en el que

ya que la ecuación anterior se convierte en

Es posible llegar a este mismo resultado analizando el diagrama de estimación por mínimos cuadrados (Figura 9.1 de Isermann) como se sigue.

y tenemos que

. Nótese que para obtener el error de predicción, el diagrama nos sugiere hacer

. Tal y como en el caso anterior,

, por tanto

Estas últimas ecuaciones tendrán repercursiones fundamentales cuando estudiemos la convergencia de los parámetros a sus valores reales.

Ahora bien, para los análisis anteriores solo hemos utilizado un único valor para la predicción. En general, necesitaremos recopilar muchos más datos (como mínimo

. Para simplificar cálculos en adelante, asumiremos que

. Podemos plantear entonces el vector de observaciones

y una matriz de predictores asociados

Para llegar al problema de regresión lineal dado por

que resumiremos de la forma

. Este problema, en general, no tendrá solución, razón por la cual recurriremos a la pseudo-inversa para llegar al problema

Este problema podemos abordarlo desde la perspectiva de las funciones de costo o directamente desde el álgebra lineal. Para el caso de la función de costo, tenemos que minimizar

Notando que

se llega finalmente a que

Para el caso de la perspectiva del álgebra lineal, las deducciones son mucho más sencillas. Solo resta pedir que

, siendo

el espacio columna de Φ. De esta forma, se sigue inmediatamente que

. Para el caso de la segunda condición basta con que

puesto que

El método de mínimos cuadrados no recursivo puede establecerse como:

Recopilar los datos del modelo (N+m+d).
Construir las matrices .
Determinar los elementos de .
Resolver .

Es particularmente interesante estudiar los elementos que componen a

y a

Nótese que

y por lo tanto

. De tal modo, podemos llamar de ahora en adelante

Siguiendo la misma lógica del empleada anteriormente, resulta fácil ver que

, y llamamos

Por tanto

Estudio de la convergencia a los parámetros reales en el método de mínimos cuadrados

Nos preguntamos ahora por la esperanza de los parámetros estimados

Recordando que

llegamos a que

Ahora, recordando que

, llegamos a que

Estudiemos la expresión

, es decir que

Si este vector es 0, entonces el sesgo desaparece y los parámetros se estiman correctamente. Para cumplir tal condición necesitaríamos que

Pero esta es la ecuación de Yule-Walker para una señal auto-regresiva y se cumple cuando

, siendo

un ruido blanco, es decir, la estimación será insesgada cuando el ruido coloreado que ingresa al sistema se produce por acción del filtro

. Es decir la estructura del modelo es ARX.

Esta consecuencia nos muestra inmediatamente una forma de validar nuestra estimación:

y también teníamos que

, por tanto

. Asumiendo que los parámetros se estiman correctamente llegamos, finalmente, a que

, es decir, el error de las predicciónes de nuestro modelo debe corresponder con un ruido blanco.

Los resultados anteriores se resumen en la siguiente gráfica:

Manejo de la incertidumbre en la estimación de los parámetros

Condiciones para la identificación del modelo

Método de mínimos cuadrados recursivo

Clase 8: Método de mínimos cuadrados (aplicaciones)

Clase 9: Métodos de variable instrumental

Clase 10: Métodos de variable instrumental

Funciones adicionales para el curso

function funresult=ldiv(a,b,N)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%Function (ldiv) : calculate inverse Z-transform by long division 
% Author : Tamer mohamed samy abdelazim Mellik
% Contact information : 
%Department of Electrical & Computer Engineering,
%University of Calgary,
%2500 University Drive N.W. ,
%Calgary, AB T2N 1N4 ,
%Canada .
% email :abdelasi@enel.ucalgary.ca  
% email : tabdelaz@ucalgary.ca
% Webpage : http://www.enel.ucalgary.ca/~abdelasi/
% Date    : 2-5-2002
% Version : 1.0.0
%Example
% This function like deconv but it help if the numerator less or equal degree of denominator
% if you have this function (It must arranged in terms of minus power of Z):
%              1
% G(z)= -----------------      
%             -1     -2
%      ( 5 - Z  - 3 Z  )
% and you want to calculate long division or inverse Z transform :
% The numerator  is a=[1] and  the denominator is b= [5 -1 -3 ]
% call the function ldiv(a,b) to get the funresult 20 items (default)
% another example :
%                -2     -3 
%       ( 5 - 3 Z  + 4 Z  )
%  G(z)-----------------      
%              -1     -2
%        (5 - Z  - 3 Z )
% a=[5 0 -3 4] , b= [5 -1 -3 ] and you want the funresult 100 terms !
% ldiv(a,b,100) 
% Note : The author doesn't have any responsibility for any harm caused by the use of this file
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%a numerator
%b denominator
%default order of the filter == 20
funresult=[];
if nargin < 3
    if nargin > 1
            N=20;
    else
            disp('Usage: M = ldiv(a,b,N)')
            disp('a:numerator , b denominator and N is the order of the resultant filter')
            return
    end 
end
if size(a) < 1
    disp('Error: numerator must at least have one element not empty')
    return
end
if size(b) < 1
    disp('Error: denominator must at least have one element not empty')
    return
end
if b(1)==0
    disp('Error: The first element of denominator must have nonzero value')
    return
end
if size(b) < 2
    funresult=a./b;
    for i =length(funresult)+1:N
       funresult(i)=0;
    end
    return
end
for i = length(a)+1:N
    a(i)=0;
end
for i = 1 : N
    funresult(i)=a(1)/b(1);
    if length(a)>1
            for k= 2:length(b)
                if k > length(a)
                    a(k)=0;
                end
                a(k)=a(k)-funresult(length(funresult))*b(k);
            end
        for i = 1:length(a)-1
                a(i)=a(i+1);
        end
        a=a(1:length(a)-1);
    end
end
end