Aproximación de una f. de transferencia por minimos cuadrados recursivos

Contents

Cáculo de la matriz de coeficientes theta

M=[0 0 0 0;1 0 0 0;1 1 -0.99 0; 1 1 -1.1 -0.99; 1 1 -0.46 -1.1;1 1 -0.66 -0.46; 1 1 -0.96 -0.66]
Y=[0;0.99;1.1;0.46;0.66;0.96;0.72]
theta_0=inv((M'*M))*M'*Y

M =

         0         0         0         0
    1.0000         0         0         0
    1.0000    1.0000   -0.9900         0
    1.0000    1.0000   -1.1000   -0.9900
    1.0000    1.0000   -0.4600   -1.1000
    1.0000    1.0000   -0.6600   -0.4600
    1.0000    1.0000   -0.9600   -0.6600


Y =

         0
    0.9900
    1.1000
    0.4600
    0.6600
    0.9600
    0.7200


theta_0 =

    0.9900
    0.5192
    0.4129
    0.5994

Cáculo de la nueva theta al ampliar valores

%Ampliación trabajando desde cero, es decir, ejecutando el proceso completo
m=[1 1 -0.72 -0.96]
Y1=0.63
Mampliada=[M;m]
Yampliada=[Y;Y1]
theta_1=inv((Mampliada'*Mampliada))*Mampliada'*Yampliada

%Ampliación mediante metodo de recursividad
PK=inv(M'*M)
error=Y1-(m*theta_0)
PK1=PK*(diag(ones(1,4))-((m'*m*PK)/(1+m*PK*m')))
theta_2=theta_0+PK1*m'*error

%Comprobación. La resta debería ser cero
fprintf('restamos theta_1 y theta_2. Si da cero es que los dos metodos dan el mismo resultado')
theta_1-theta_2

m =

    1.0000    1.0000   -0.7200   -0.9600


Y1 =

    0.6300


Mampliada =

         0         0         0         0
    1.0000         0         0         0
    1.0000    1.0000   -0.9900         0
    1.0000    1.0000   -1.1000   -0.9900
    1.0000    1.0000   -0.4600   -1.1000
    1.0000    1.0000   -0.6600   -0.4600
    1.0000    1.0000   -0.9600   -0.6600
    1.0000    1.0000   -0.7200   -0.9600


Yampliada =

         0
    0.9900
    1.1000
    0.4600
    0.6600
    0.9600
    0.7200
    0.6300


theta_1 =

    0.9900
    0.5185
    0.4119
    0.6012


PK =

    1.0000   -1.0000    0.0000    0.0000
   -1.0000    5.3118    3.7543    1.5276
    0.0000    3.7543    3.9360    0.7347
    0.0000    1.5276    0.7347    1.4250


error =

   -0.0065


PK1 =

    1.0000   -1.0000    0.0000    0.0000
   -1.0000    5.2967    3.7315    1.5667
    0.0000    3.7315    3.9015    0.7939
    0.0000    1.5667    0.7939    1.3233


theta_2 =

    0.9900
    0.5185
    0.4119
    0.6012

restamos theta_1 y theta_2. Si da cero es que los dos metodos dan el mismo resultado
ans =

  1.0e-014 *

    0.0444
   -0.1887
   -0.0278
   -0.1221

Cálculo de la función de coste

%Creamos los vectores uk_2 (entradas empezando en k-2) e yk_2 (salidas empezando en k-2)
uk_2=[0 0 1 1 1 1 1 1];
yk_2=[0 0 0 0.99 1.1 0.46 0.66 0.96];

%Calculamos la función de coste del primer caso
coste=0;
j=3;
for i=1:1:length(Y)
    coste=coste+(Y(i)-(uk_2(j-1)+theta_0(2)*uk_2(j-2)-theta_0(3)*yk_2(j-1)-theta_0(4)*yk_2(j-2))^2);
    j=j+1;
end
fprintf('el coste en el primer caso es:')
disp(coste)

%Calculamos la función de coste con la ampliación
uk_2=[0 0 1 1 1 1 1 1 1];
yk_2=[0 0 0 0.99 1.1 0.46 0.66 0.96 0.72];

coste=0;
j=3;
for i=1:1:length(Yampliada)
    coste=coste+(Yampliada(i)-(uk_2(j-1)+theta_1(2)*uk_2(j-2)-theta_1(3)*yk_2(j-1)-theta_1(4)*yk_2(j-2))^2);
    j=j+1;
end
fprintf('el coste en el segundo caso es:')
disp(coste)

el coste en el primer caso es:    0.5143

el coste en el segundo caso es:    0.7347

Representación gráfica

%Creamos el vector de tiempos y representamos los datos del ejercicio en
%azul
k=[0 1 2 3 4 5 6];
figure(1)
hold on
plot(k,Y)%en azul los datos del ejercicio

%Representamos en rojo los resultados generados por la función de
%transferencia obtenida en la primera aproximación
num = [1 theta_0(2,1)];
den = [1 theta_0(3,1) theta_0(4,1)];
G_0 = filt(num,den)
step(G_0,'red')

%Representamos en verde los resultados generados por la función de
%transferencia obtenida en la primera aproximación
num = [1 theta_1(2,1)];
den = [1 theta_1(3,1) theta_1(4,1)];
G_1 = filt(num,den)
step(G_1,'green')

hold off
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

 
Transfer function:
       1 + 0.5192 z^-1
-----------------------------
1 + 0.4129 z^-1 + 0.5994 z^-2
 
Sampling time: unspecified
 
Transfer function:
       1 + 0.5185 z^-1
-----------------------------
1 + 0.4119 z^-1 + 0.6012 z^-2
 
Sampling time: unspecified

Published with MATLAB® 7.10