• Author / Uploaded
• dante javier solis oncoy

Citation preview

Cuarta Practica Calificada Domiciliaria En la figura se muestra un sistema de control de nivel del líquido h(t) de un recipiente cilíndrico de área Ah, altura H y resistencia hidráulica R h en el caudal de salida. El caudal de entrada q e(t) es gobernado por una electroválvula que consiste en un impulsor hidráulico (de desplazamiento x(t), masa M, fricción B, resorte de constante K y sección de presión A p) y un conversor voltaje / presión. El sensor de nivel consiste en flotador el cual está unido a una cremallera de carrera vertical que al desplazarse hace girar a un potenciómetro rotatorio de 3 vuelta y radio R.

Figura 01 Sistema de Control de nivel de un tanque

Con las ecuaciones físicas adicionales del sistema:

p ( t ) =Cec ( t ); f ( t )= A p p ( t ); τ q q˙ e ( t ) +q e ( t ) =k q x (t ); e h (t )=k ϑ θ (t) Datos de los parámetros: H=2m; Ah=4m2; Rh =60seg/m2; Ap.=0.0025m2; M = 0.25kg; B = 0.075 N seg/m; K = 0.2N/m; C=8V/m 2; kq= 0.05m2; τq =2seg; R = 0.10611m; kθ = 5/3π vol/rad

Para el Nro. de Orden igual a 28, se tiene:

t s=

N 0 +19 =23.5 minutos=1410 segundos 2

T=

0.5 N 0+ 6 =1.25 16

Se tendrá que utilizar un Integrador en atraso y el método de discretización de la planta será tipo “tustin”.

Figura 02 Sistema de control de seguimiento con realimentación de estado y control integral. K=matriz de ganancia de realimentación de estado Ki=ganancia integral del error

Figura 03 Sistema de control de seguimiento con realimentación de estado estimado y control integral K=matriz de ganancia del controlador Ke=matriz de ganancia del estimador Ki=ganancia integral del error

Usando el diagrama de la figura (01) y su respectiva representación en el espacio de estado en las partes (b) y (c) de la tercera práctica se pide: a) Realizar la conversión de tiempo continuo a tiempo discreto del espacio de estado

representado en la parte (b) y representarlo en un diagrama de bloques de estado en tiempo discreto con un periodo de muestreo T si se colocan muestreadores sincronizados δ m (t) en e h (t) y en e c (t) y un retenedor de orden cero después del muestreador en e c (t). Vamos a considerar el gráfico anterior como:

Agregando los muestreadores y el retenedor de orden cero:

Entonces en tiempo discreto usando el espacio de estados llegaremos a:

x ( k +1 )=Gx ( k )+ Hu ( k ) … .(1) y ( k )=Cx ( k ) + Du(k ) De la practica 3 se obtuvo la ecuación en espacio de estados:

x˙ 1 ( t ) −0.00416 0.25 0 0 x˙ 2 ( t ) = 0 −0.5 0.025 0 0 0 0 1 x˙ 3 ( t ) 0 0 −0.8 −0.3 x˙ 4 (t )

[ ][

x1 ( t ) 0 x 2 ( t ) + 0 u(t) 0 x3 ( t ) 0.08 x 4 (t )

][ ] [ ]

x1 ( t ) y ( t ) =[ 1 0 0 0 ] x 2 ( t ) + [ 0 ] u(t ) x3 ( t ) x4 ( t )

[]

Siendo las siguientes matrices:

−0.00416 0.25 0 0 0 −0.5 0.025 0 , B= A= 0 0 0 1 0 0 −0.8 −0.3

[

0 0 , C=[ 1 0 0 0 ] y D= [ 0 ] 0 0.08

][]

Utilizando el siguiente comando en MATLAB para el método Tustin: clc clear all close all %discretizacion en espacio de estados(tustin) A=[-0.00416 0.25 0 0;0 -0.5 0.025 0;0 0 0 1; 0 0 -0.8 -0.3]; B=[0;0;0;0.08]; C=[1 0 0 0]; D=[0]; Ts=1.25; [G,H,C1,D1]=c2dm(A,B,C,D,Ts,'tustin')

Tenemos las siguientes matrices:

0.9948 0.2375 0.0029 0.0015 0 0.5238 0.01880.0099 , G= 0 0 0.5833 0.8333 0 0 −0.6667 0.3333

[

]

0.0001 0.0005 H= 0.0417 0.0667

[ ]

C1=[ 0.9974

0.1187 0.0015 0.0008 ]

D1=[ 3.8652e-05 ]= [ 0.000038652 ]

El diagrama de bloques será:

b) Determinar la matriz de transferencia pulso en lazo abierto del sistema para este

periodo T Tomando transformada Z a la ecuación (1):

zX ( z )−zx ( 0 )=GX ( z ) + HU ( z )

Y ( z )=CX ( z ) + DU (z) Si consideramos las condiciones iníciales nulas:

X ( z )=( zI−G )−1 HU ( z) Y ( z )=[ C ( zI−G )−1 H + D ] U (z) Luego, La matriz de transferencia vendrá dada por:

F ( z )=C ( zI −G )−1 H + D …( 19) Dada la similitud que hay com esta matriz de transferencia en tiempo discreto con la de tiempo continuo, podemos usar el mismo comando para el cálculo de esta matriz: Utilizando el siguiente comando en MATLAB: clc clear all close all A=[-0.00416 0.25 0 0;0 -0.5 0.025 0;0 0 0 1; 0 0 -0.8 -0.3]; B=[0;0;0;0.08]; C=[1 0 0 0]; D=[0]; Ts=1.25; [G,H,C,D]=c2dm(A,B,C,D,Ts,'tustin'); [NUZ,DEZ] = ss2tf(G,H,C,D); tf(NUZ,DEZ); R=tf(NUZ,DEZ,Ts) Se tiene la siguiente función de transferencia: 4 3 2 Y (z) −3 0.03865 z +0.1546 z + 0.2319 z +0.1546 z+ 0.03865 =10 U ( z) z 4−2.435 z 3+ 2.663 z 2−1.617 z +0.3908

c) Diseñar un sistema de control de seguimiento con realimentación de estado de ganancia K y control integral de ganancia K i usando Control Óptimo Cuadrático en tiempo discreto según figura 02, eligiendo adecuadamente las matrices de ponderación Q y R para que no haya desbordamiento en el tanque. Se tiene el diagrama:

De la figura se obtienen las siguientes ecuaciones:

x ( k +1 )=Gx ( k )+ Hu ( k ) y ( k )=Cx ( k ) + Du(k ) u ( k ) =−Kx ( k ) + K I v (k ) v ( k +1 )=v ( k ) +r (k )−K M y ( k ) Se obtiene un sistema de control ampliado dado por:

(

x ( k +1 ) G 0 x (k ) + H = u (k )+ 0 r ( k ) −K M C 1 v ( k ) −K M D 1 v ( k +1 )

)(

)( ) (

y ( k )= ( C

0)

) ()

x(k) + Du(k ) v (k )

( )

u ( k ) =−( K −K I )

( vx (( kk )))

Ahora las nuevas matrices son:

^x ( k )=

( xv (( kk )))

^= G

(−KG C 01)

^ (C 0) C= Se define:

^x e ( k )=^x ( k )−^x ( ∞ )

M

^ K= ( K −K I ) ue ( k )=u ( k )−u ( ∞ )

Al restar las ecuaciones se obtiene:

^ x^ e ( k ) + ^ ^ x^ e ( k ) ^x e ( k +1 )=G H u e ( k )ue ( k )=− K Finalmente:

^ H ^^ ^x e ( k +1 )=( G− K ) ^x e ( k ) Con ecuación característica:

^ H=

^ H ^^ |zI− G+ K|=0

(−KH D ) H =( 01) ¿

M

^ D=D r ( k +1 ) =r ( k )=r

Donde, al elegir la entrada de referencia 1 voltio r(t)=1:

e (t )=r ( t )−K M y ( t ) → e ( ∞ ) =0=r ( ∞ ) −K M y ( ∞ ) → K M =

r (∞) 1 = =0.5 y (∞ ) 2

En el problema de control óptimo cuadrático se desea determinar una ley para el vector de control u(k) tal que un índice de desempeño cuadrático se minimice. El índice de desempeño cuadrático con el que trabajaremos será: ∞

J=

1 ∑ [ ^x T ( k ) Q^ ^x e( k )+u eT (k ) R^ u e( k) ] 2 k =0 e

^ , es simétrica definida positiva o semidefinida positiva de Donde la matriz de ponderación Q 5x5 y la matriz de ponderación ^ R , es simétrica y siempre definida positiva de 1x1. La matriz de ganancia óptima en estado estacionario ^ K se determina por la siguiente expresión:

^ ^+ ^ ^ ]−1 ^ ^ … (2) K=[R HT ^ PH HT ^ PG Donde la matriz ^ P se obtiene de la ecuación de Ricatti: ^ ^ G ^T P ^ G− ^ G ^T P ^^ ^T P ^G ^ …(3) P=Q+ H[^ R+ ^ HT ^ P^ H ]−1 H La solución de la ecuación de Ricatti se obtiene mediante la ecuación: −1 ^ ^+G ^T P ^ (k ) G− ^ G ^T ^ ^[^ ^T P ^ (k ) ^ ^ … (4) P (k + 1)= Q P (k ) H R+ H H] ^ HT ^ P (k ) G

Eligiendo valores de Q y R que cumplan la condición T S=23.5min=1410 seg :

10000 0 0 0 0 0 10000 0 0 0 Q= 0 0 100 0 0 y R=( 100 ) 0 0 0 100 0 0 0 0 0 2

(

)

Calculamos las matrices P y K con el siguiente código Matlab: clc clear all close all CC=[0.9974 0.1187 0.0015 0.0008 0]; D=[0]; GG=[0.9948 0.2375 0.0029 0.0015 0;0 0.5238 0.0188 0.0099 0;0 0 0.5833 0.8333 0;0 0 -0.6667 0.3333 0;-0.4987 -0.05935 -0.00075 -0.0004 1]; HH=[0.0001;0.0005;0.0417;0.0667;-1.9325e-05]; Q=[10000 0 0 0 0;0 10000 0 0 0;0 0 10000 0 0;0 0 0 10000 0;0 0 0 0 2]; R=[100]; P=diag(0,4); for n=1:200; P=Q+GG'*P*GG-GG'*P*HH*inv(R+HH'*P*HH)*HH'*P*GG

Las matrices P Y K serán:

P= 1.0e+05 * 8.4307 4.1096 0.1260 0.0992 -0.0878 4.1096 2.1742 0.0642 0.0503 -0.0429 0.1260 0.0642 0.2174 0.0597 -0.0013 0.0992 0.0503 0.0597 0.2068 -0.0010 -0.0878 -0.0429 -0.0013 -0.0010 0.0035 Ks = 5.5948 2.8270 -1.1586 6.2300 -0.0583 Entonces: K= [5.5948 2.8270 -1.1586 6.2300 -0.0583] Ki= [-0.0583] d) Si e r ( t ) =r (k ) determinar la Matriz de Transferencia Pulso en lazo cerrado: Y(z)/R(z). De las ecuaciones anteriores:

^ ^x ( k )+ H ^ u^ ( k )+ 0 r ( k ) ^x ( k )=G 1 ^ ^x ( k )+ D ^ u^ ( k ) y ( k )= C u^ ( k ) =− ^ K ^x ( k )

()

( 01)r ( k )

^ ^ ^ ) ^x ( k ) + ^x ( k +1 )=( G− HK Usando Transformada Z:

^ ^ ^ ) ^x ( k ) y ( k )= ( C− DK

^ ^x ( z )+ ^ ^ ^x ( z )= 0 r ( z ) → ( zI −G ^ +^ ^ ) x^ ( z ) = 0 r ( z ) … (1) zI ^x ( z )−G HK HK 1 1 ^ ^x ( k )+ D ^ ^ ^ u ( k ) =( C− y ( k )= C D^ K ) x^ ( k ) …(2)

()

De las ecuaciones (1) y (2):

()

−1

^ −^ ^ +^ Y ( z ) =( C D^ K ) ( zI−G H^ K)

(01 ) R ( z )

Y ( z) ^ ^ ^ −1 ^ H ^^ =( C− D K ) ( zI−G+ K) 0 1 R( z)

()

Usamos el siguiente código de Matlab para calcular la función de transferencia: clc clear all close all CC=[0.9974 0.1187 0.0015 0.0008 0]; D=[0]; GG=[0.9948 0.2375 0.0029 0.0015 0;0 0.5238 0.0188 0.0099 0;0 0 0.5833 0.8333 0;0 0 -0.6667 0.3333 0;-0.4987 -0.05935 -0.00075 -0.0004 1]; HH=[0.0001;0.0005;0.0417;0.0667;-1.9325e-05]; Ks =[4.7955 2.4348 -1.1706 6.2205 -0.0317]; P=CC-D*Ks; Q=GG-HH*Ks; R=[0 ;0 ;0 ;0 ;1]; S=0; [NUMZ,DENZ]=ss2tf(Q,R,P,S); Ts=1.25; L=tf(NUMZ,DENZ,Ts)

Y ( z) 10−5 × ( 0.8717 z 3 +0.2736 z 2+ 0.7949 z +0.04407 ) = 5 R ( z) z −3.067 z 4 +3.713 z3 −2.367 z2 +0.8697 z−0.1483 e) Graficar las variables de estado x(k), del error e(k), v(k), y(k) y la variable de control u(k) ante una entrada en r(k) (tipo escalón) para llenarlo completamente al tanque sin desbordamiento usando como vector de condiciones iniciales x ( 0 ) el mismo de la parte (h) de la tercera práctica.

Si tomamos los valores iniciales como en la tercera práctica, X (0) =0.1 cero. Obtenemos las gráficas:

m3 , y llenando el tanque desde s

300

250

200

150

100

50

0

-50 0

500

1000

1500

2000

2500

3000

3500

4000

4500

RESPUESTA X1(t) [RESPUESTA h]

2.5

ALTURA (metros)

2

1.5

1

0.5

0 0

500

1000

1500

2000

2500

3000

3500

4000

4500

3500

4000

4500

TIEMPO (segundos)

RESPUESTA X2(t) [RESPUESTA qe]

0.1 0.09

0.07 0.06

3

/segundos)

0.08

CAUDAL (metros

0.05 0.04 0.03 0.02 0.01 0 0

500

1000

1500

2000

2500

TIEMPO (segundos)

3000

RESPUESTA X3(t) [RESPUESTA x]

0.7

0.6

POSICION (metros)

0.5

0.4

0.3

0.2

0.1

0

-0.1 0

500

1000

1500

2000

2500

3000

3500

4000

4500

3500

4000

4500

3500

4000

4500

TIEMPO (segundos)

RESPUESTA X4(t) [RESPUESTA dx/dt]

0.01

0

-0.005

-0.01

-0.015

-0.02 0

500

1000

1500

2000

2500

3000

TIEMPO (segundos)

RESPUESTA X5(t) [RESPUESTA v]

300

250

ENTRADA PROPORCIONAL (metros)

VELOCIDAD (metros/segundos)

0.005

200

150

100

50

0 0

500

1000

1500

2000

2500

TIEMPO (segundos)

3000

GRAFICA DE e(kT)

1 0.8 0.6 0.4

e(kT)

0.2 0 -0.2 -0.4 -0.6 -0.8 -1 0

500

1000

1500

2000

2500

3000

3500

4000

4500

3000

3500

4000

4500

3000

3500

4000

4500

TIEMPO (segundos)

GRAFICA DE y(kT)

4

3.5

3

y(kT)

2.5

2

1.5

1

0.5

0 0

500

1000

1500

2000

2500

TIEMPO

GRAFICA U(kT)

7

6

5

u(kT)

4

3

2

1

0 0

500

1000

1500

2000

2500

TIEMPO (segundos)

Para las siguientes gráficas se uso el siguiente código de Matlab:

clc close all Q=[10000 0 0 0 0;0 10000 0 0 0;0 0 10000 0 0;0 0 0 10000 0;0 0 0 0 1]; R=[100]; CC=[0.9974 0.1187 0.0015 0.0008 0]; D=[0]; GG=[0.9948 0.2375 0.0029 0.0015 0;0 0.5238 0.0188 0.0099 0;0 0 0.5833 0.8333 0;0 0 -0.6667 0.3333 0;-0.4987 -0.05935 -0.00075 -0.0004 1]; HH=[0.0001;0.0005;0.0417;0.0667;-1.9325e-05]; P=diag(0,4); for n=1:200; P=Q+GG'*P*GG-GG'*P*HH*inv(R+HH'*P*HH)*HH'*P*GG; end; Ks=inv(R+HH'*P*HH)*HH'*P*GG; M=CC-D*Ks; N=GG-HH*Ks; R=[0 ;0 ;0 ;0 ;1]; S=0; [NUMZ,DENZ]=ss2tf(N,R,M,S); Ts=1.25; L=tf(NUMZ,DENZ,Ts); X0=[0 0.1 0 0 0]; U=ones(3500,1); [Y,X]=dlsim(N,R,M,S,U,X0); stairs((1:3500)*1.25,X);grid on; figure(2);stairs((1:3500)*1.25,X(:,1));grid on; title('RESPUESTA X1(t) [RESPUESTA h]'); xlabel('TIEMPO (segundos)'); ylabel('ALTURA (metros)'); figure(3);stairs((1:3500)*1.25,X(:,2));grid on; title('RESPUESTA X2(t) [RESPUESTA qe]'); xlabel('TIEMPO (segundos)'); ylabel('CAUDAL (metros^3/segundos)'); figure(4);stairs((1:3500)*1.25,X(:,3));grid on; title('RESPUESTA X3(t) [RESPUESTA x]') ;xlabel('TIEMPO (segundos)'); ylabel('POSICION (metros)'); figure(5);stairs((1:3500)*1.25,X(:,4));grid on; title('RESPUESTA X4(t) [RESPUESTA dx/dt]'); xlabel('TIEMPO (segundos)'); ylabel('VELOCIDAD (metros/segundos)'); figure(6);stairs((1:3500)*1.25,X(:,5));grid on; title('RESPUESTA X5(t) [RESPUESTA v]'); xlabel('TIEMPO (segundos)'); ylabel('ENTRADA PROPORCIONAL (metros)'); figure(7);r=ones(3500,1);E=r-Y;stairs((1:3500)*1.25,E);grid on; title('GRAFICA DE e(kT)'); xlabel('TIEMPO (segundos)');

f)

Interpretar físicamente la dinámica del sistema basándose en las gráficas de la parte (e).

Se puede observar de las gráficas que para una condición inicial de 0.1 para el caudal de entrada que tomamos para la práctica 3; el caudal de entrada, la altura de agua del tanque se eleva al comienzo pero debido al control realizado en el tanque, éste se llena igualmente para el valor requerido de aproximadamente 1400 segundos, además que el error en estado estacionario se hace cero en el infinito, el tanque se llena de agua a la altura de 2 metros sin que exista desbordamiento. Análogamente las gráficas de los otros elementos como el caudal, la posición y la velocidad de la válvula que controla la salida de caudal trabajan de forma correcta ya que no se ve ningún sobrepico; a la vez que sus curvas son suaves y todas llegan a establecerse para un determinado tiempo; para nuestro caso aprox. a los 1400 segundos. g) Diseñar un estimador de estados de orden completo para ser utilizados como realimentación de estado en el sistema de control de la parte (c) según la figura 03,

x (k ) tienda a cero a una velocidad igual a 3 de tal manera que el error e E ( k )=x ( k )−~ veces o mayor que la velocidad de respuesta de la salida del sistema, usando el método de la ubicación de polos para el diseño del estimador Ke. De la figura 03:

Cuando el vector real x ( k ) no está disponible o no es medible en forma directa se diseña un estimador de estados para un vector estimado ~ x (k ) sea lo mas cercano al vector de estados real x ( k ) . Sea entonces el sistema de control de simple entrada descrito por las ecuaciones:

x ( k +1 )=Gx ( k )+ Hu ( k ) y ( k )=Cx ( k ) + Du(k ) Y teniendo en cuenta que las señales y ( k ) y ~ y (k ) se pueden medir, por lo tanto la ecuación dinámica del observador de estados, se escribe como:

Donde Ke es una matriz de ponderación de estados de dimensión nx1. Del método algorítmico tenemos que:

Donde Q=WN Además:

Calculando N y W con el siguiente código de Matlab: clc clear all close all A=[-0.00416 0.25 0 0;0 -0.5 0.025 0;0 0 0 1; 0 0 -0.8 -0.3]; B=[0;0;0;0.08]; C=[1 0 0 0]; D=[0]; Ts=1.25; [G,H,C,D]=c2dm(A,B,C,D,Ts,'tustin'); N=[C ; C*G ; C*G^2 ; C*G^3] syms z; I=round (G*G^(-1)); vpa(det(z*I - G),6)

N= 1.0000

0

0

0

0.9948 0.2317 0.0036 0.0015 0.9897 0.3545 0.0087 0.0078 0.9845 0.4191 0.0092 0.0152 ans = z^4 - 2.4359*z^3 + 2.66316*z^2 - 1.61664*z + 0.39082 Entonces:

a 1=−2. 4359 , a2=2. 66316 , a3=−1.61664 ,a 4=0.3 9082

De la definición de W:

−1. 61664 2. 66316 −2. 4359 1 2.66316 −2. 4359 10 W= −2. 4359 1 0 0 1 0 0 0

[

]

Con el siguiente código de Matlab calculamos las matrices restantes: clc clear all close all A=[-0.00416 0.25 0 0;0 -0.5 0.025 0;0 0 0 1; 0 0 -0.8 -0.3]; B=[0;0;0;0.08]; C=[1 0 0 0]; D=[0]; Ts=1.25; [G,H,C,D]=c2dm(A,B,C,D,Ts,'tustin'); N=[C ; C*G ; C*G^2 ; C*G^3] W =[-1.45 2.3655 -2.27891 1;2.3655 -2.27891 1 0;-2.27891 1 0 0;1 0 0 0] Q= W*N Q1=Q' M = [H G*H (G^2)*H (G^3)*H] T = M*W

N= 0.9974 0.9922 0.9871 0.9820

0.1187 0.2991 0.3923 0.4399

0.0015 0.0055 0.0090 0.0085

0.0008 0.0042 0.0105 0.0164

W= -1.6166 2.6632 -2.4359 1.0000 2.6632 -2.4359 1.0000 0 -2.4359 1.0000 0 0 1.0000 0 0 0

Q= -0.3924 1.2264 -1.4373 0.9974

0.0888 -0.0200 0.0098 0.1187

-0.0010 -0.0005 0.0019 0.0015

0.0008 0.0023 0.0023 0.0008

Q1 = -0.3924 0.0888 -0.0010 0.0008

1.2264 -0.0200 -0.0005 0.0023

-1.4373 0.0098 0.0019 0.0023

0.9974 0.1187 0.0015 0.0008

M= 0.0001 0.0005 0.0417 0.0667

0.0004 0.0017 0.0799 -0.0056

0.0010 0.0023 0.0420 -0.0551

0.0016 0.0015 -0.0214 -0.0463

T= 0.0001 -0.0005 0.0217 -0.0347

0.0002 -0.0005 -0.0416 0.1360

0.0002 0.0005 -0.0216 -0.1679

0.0001 0.0005 0.0417 0.0667

Calculando los valores de ∝n: Calculamos los polos deseados, que aparentemente se usaron en el control óptimo usando el método de ubicación de polos, para poder escoger de forma arbitraria los polos deseados para el estimador tal que la velocidad de respuesta de la salida sea 3 veces más que para el control óptimo. Para lograr esto se tomará la matriz de ganancia K hallada en la parte (c) y se hallarán los polos deseados:

K= [α4-a4

α3-a3

α2-a2

α1-a1]T −1

Entonces:

[ α4-a4

α3-a3

Usando Matlab calculamos Z=KT: z= -0.2423 0.8953 -1.0185 0.3689

α2-a2

α1-a1]=KT

Entonces:

∝1=−2.067 , ∝2=1.6446 , ∝3 =−0.72134 , ∝4=0.14852 Teniendo las raíces de las condiciones iniciales tenemos que tratar de triplicar la velocidad de respuesta de la salida. Usamos el siguiente comando de Matlab para graficar el error:

clc clear all Q=[10000 0 0 0 0;0 10000 0 0 0;0 0 10000 0 0;0 0 0 10000 0;0 0 0 0 2]; R=[100]; CC=[0.9974 0.1187 0.0015 0.0008 0]; D=[0]; GG=[0.9948 0.2375 0.0029 0.0015 0;0 0.5238 0.0188 0.0099 0;0 0 0.5833 0.8333 0;0 0 -0.6667 0.3333 0;-0.4987 -0.05935 -0.00075 -0.0004 1]; HH=[0.0001;0.0005;0.0417;0.0667;-1.9325e-05]; P=diag(0,4); for n=1:200; P=Q+GG'*P*GG-GG'*P*HH*inv(R+HH'*P*HH)*HH'*P*GG; end; Ks=inv(R+HH'*P*HH)*HH'*P*GG; M=CC-D*Ks;N=GG-HH*Ks; R=[0 ;0 ;0 ;0 ;1];S=0;[NUMZ,DENZ]=ss2tf(N,R,M,S); Ts=0.617647;L=tf(NUMZ,DENZ,Ts); step(L) grid on title('RESPUESTA A UN ESCALON PARA NUESTRO SISTEMA') xlabel('TIEMPO (segundos)') ylabel('ALTURA (metros)') A=GG'-CC'*Ks; B=CC'; C=HH'; D=0; [NUZ,DEZ] = ss2tf(A,B,C,D); T=0.617647; N2=tf(NUZ,DEZ,T); step(N2) grid on title('RESPUESTA A UN ESCALON PARA NUESTRO SISTEMA') xlabel('TIEMPO (segundos)') ylabel('ALTURA (metros)')

h) Graficar las variables de estado x(k), ~ x (k ), del error e(k), v(k), y(k) y la variable de control u(k) ante una entrada en r(k)=1 rad (tipo escalón) usando como vector de condiciones iniciales e E ( 0 )=( 0.5 i) j)

0.2 −0.01 0 )T .

Comparar e interpretar físicamente la respuesta del sistema basándose en las gráficas de la parte (e) y (h). Determinar la Función de Transferencia Pulso del controlador-estimador U(z)/-Y(z).