Tarea 3: Diseño de Controladores y Observadores de Estado para un Sistema MIMO Sea el siguiente proceso multivariable multivariable descrito por las ecuaciones diferenciales:
Donde
,
⃛ ̈ ̈ ̇ ⃛ ̈ ̈ ̇ ̇ ̈ ̇ ̈
: las variables de salida;
: las variables de entrada
a) Eligiendo como variable de estado
Representaremos al sistema mediante variables de estado en tiempo continuo y con un muestreador de periodo T=0.04 seg y un retenedor de orden cero, discretizaremos el proceso. Llevando a matrices el siguiente sistema, tendremos:
̇̇ ̇̇ [̇̇ ] [ ] [] [ ] [] ̇
En forma compacta el proceso se puede representar como
El tiempo discreto con muestreador de periodo T=0.04 seg. Las ecuaciones en el espacio de estado se convierten en
Usaremos el MATLAB para convertir el sistema a Digital y hallar las matrices G y H.
Programa en MATLAB % a) Sistema continuo, luego retenedor de Orden Zero A=[0 1 0 0 0 0;0 0 1 0 0 0;-8 -3 -4 -6 0 -5;0 0 0 0 1 0;0 0 0 0 0 1;-7 0 -6 -5 -2 4] B=[0 0 0;0 0 0;8 6 -2;0 0 0;0 0 0;6 8 -3] C=[1 0 0 0 0 0;0 0 0 1 0 0] D=[0 0 0;0 0 0] I=[1 0 0 0 0 0;0 1 0 0 0 0;0 0 1 0 0 0;0 0 0 1 0 0;0 0 0 0 1 0;0 0 0 0 0 1] Ts=0.04 [Gd Hd Cd Dd]=c2dm(A,B,C,D,Ts,'zoh' Dd]=c2dm(A,B,C,D,Ts,'zoh') )
Respuesta del MATLAB >> Planta_SS Gd = 0.9999 0.0400 0.0008 -0.0001 0.0000 -0.0001 -0.0057
0.9976
0.0372 -0.0043
0.0000 -0.0040
-0.2698 -0.1174
0.8728 -0.2033
0.0037 -0.2023
-0.0001 0.0000 -0.0001 0.9999 0.0400 0.0008 -0.0054 0.0001 -0.0048 -0.0039 0.9983 0.0437 -0.2672
0.0091 -0.2427 -0.1895 -0.0912
Hd = 0.0001
0.0001 -0.0000
0.0058
0.0041 -0.0014
0.2738
0.1913 -0.0624
0.0001
0.0001 -0.0000
0.0046
0.0064 -0.0024
0.2235
0.3205 -0.1214
Cd = 1
0
0
0
0
0
0
0
0
1
0
0
1.1971
Dd = 0
0
0
0
0
0
De aquí obtenemos las matrices Gd, Hd.
[ [
] ]
b) Mediante el método de ubicación de polos para un sistema multivariable, diseñaremos un controlador en realimentación de estados de ganancia K para que el vector de salida y(kT) siga un vector de referencia r(kT) en estado estacionario corrigiendo el error en estado estacionario mediante una matriz de prealimentación
, con la siguiente ley de control
. Las variables de respuesta deben ser moderadas en el tiempo
transitorio y apropiadas en régimen permanente para que las variables de control u(kT) no presenten grandes magnitudes. Vamos a ver si el Sistema es Completamente controlable, esto significa que la matriz de controlabilidad M deberá tener un rango=6.
Ingresando la ecuación para hallar M en MATLAB tenemos: M=[Hd Gd*Hd Gd^2*Hd Gd^3*Hd Gd^4*Hd Gd^5*Hd] rango_M=rank(M)
El MATLAB nos devuelve: M= Columns 1 through 7
0.0001 0.0001 -0.0000 0.0005 0.0004 -0.0001 0.0012 0.0058 0.0041 -0.0014 0.0151 0.0100 -0.0032 0.0214 0.2738 0.1913 -0.0624 0.1931 0.1016 -0.0298 0.1260 0.0001 0.0001 -0.0000 0.0004 0.0006 -0.0002 0.0011 0.0046 0.0064 -0.0024 0.0130 0.0195 -0.0074 0.0208 0.2235 0.3205 -0.1214 0.2007 0.3367 -0.1300 0.1922
Columns 8 through 14 0.0008 -0.0003 0.0022 0.0013 -0.0004 0.0032 0.0017 0.0124 -0.0038 0.0253 0.0115 -0.0031 0.0269 0.0074 0.0193 0.0007 0.0681 -0.0612 0.0311 0.0157 -0.1452 0.0017 -0.0006 0.0021 0.0033 -0.0013 0.0034 0.0057 0.0336 -0.0129 0.0286 0.0499 -0.0194 0.0368 0.0695 0.3766 -0.1477 0.1972 0.4426 -0.1756 0.2162 0.5393
Columns 15 through 18 -0.0005 0.0043 0.0018 -0.0005 -0.0013 0.0265 -0.0002 0.0020 0.0633 -0.0346 -0.2380
0.0995
-0.0022 0.0050 0.0089 -0.0035 -0.0271 0.0461 0.0936 -0.0368 -0.2157 0.2505 0.6731 -0.2706
rango_M = 6
Ahora usaremos la matriz aleatoria Ks:
La Matriz
Usando el siguiente código en Matlab: Ks=[1 -3 4 5 6 -5;4 13 -4 2 -2 4;1 -2 4 0 -3 2] G0=Gd-Hd*Ks m=eig(G0)
Matlab nos devuelve: Ks = 1
-3
4
4
13
-4
1
-2
4
5 2 0
6
-5
-2 -3
4 2
G0 = 0.9996 0.0394 0.0008 -0.0006 -0.0004 0.0001 -0.0267
0.9584
0.0362 -0.0415 -0.0304
0.0110
-1.2463 -1.9074
0.7924 -1.9549 -1.4438
0.5265
-0.0004 -0.0010 0.0002 0.9995 0.0397 0.0009 -0.0332 -0.0741 0.0121 -0.0395 0.9764 0.0458 -1.6515 -3.7303 m= 0.3798 1.2977 + 0.2744i 1.2977 - 0.2744i 0.9625 0.9925 1.0713
0.6311 -1.9481 -1.1555
1.2754
Vemos que los valores propios de la matriz
Sea
son distintos, esto es:
entonces la matriz
El código en Matlab sería: v=[1 -4 2]' H0=Hd*v
Matlab nos regresa: v= 1 -4 2 H0 = -0.0002 -0.0135 -0.6161 -0.0003 -0.0258 -1.3015 Sean los polos en tiempo continuo:
En tiempo discreto los polos serán (usando Matlab): s=[-5+5i -5-5i -3+3i -3-3i -10 -12] z0=exp(Ts*s)
Matlab nos regresa: s= -5.0000 + 5.0000i -5.0000 - 5.0000i -3.0000 + 3.0000i -3.0000 - 3.0000i -10.0000
-12.0000
z0 = 0.8024 + 0.1627i 0.8024 - 0.1627i 0.8805 + 0.1062i 0.8805 - 0.1062i 0.6703 Las matrices
representan un sistema de simple entrada y con el espectro de polos en lazo
cerrado se obtiene la matriz Luego de hallar
0.6188
.
hallaremos la matriz deseada de realimentación
Código en MATLAB: K0=place(G0,H0,z0) K=Ks+v*K0 eig(Gd-Hd*K)
Matlab nos devuelve: K0 = 37.9981
6.7827 -0.3967 -52.6251 -6.4576 -0.7807
K= 38.9981
3.7827
3.6033 -47.6251 -0.4576 -5.7807
-147.9924 -14.1310 -2.4133 212.5003 23.8306 76.9962 11.5655
3.2066 -105.2501 -15.9153
7.1229 0.4385
ans = 0.6188 0.6703 0.8024 + 0.1627i 0.8024 - 0.1627i 0.8805 + 0.1062i 0.8805 - 0.1062i Entonces la matriz K:
Vemos que los autovalores coinciden con los polos de lazo cerrado propuestos.
La matriz de Prealimentación
()
tiene dimensión 3x2 se determina utilizando la siguiente
representación obtenida en estado estacionario
Con la matriz D=0 se consigue:
La matriz
de dimensión 2x3, por lo tanto no es posible despejar la matriz
directamente, entonces se debe determinar la pseudoinversa de la matriz F que viene a ser:
Evaluando en Matlab se obtiene: F=Cd*(I-Gd+Hd*K)^-1*Hd Kr=F'*(F*F')^-1
Matlab nos devuelve: F= 0.1414
0.0299 -0.0027
0.0991
0.0248 -0.0036
Kr = 40.4359 -47.6288 -151.6189 216.5228 67.8678 -96.1976 Entonces tendríamos las matrices F y Kr:
Tambien es posible determinar la matriz de realimentación de estados utilizando la función “place” del MATLAB directamente: K=place(Gd,Hd,z0) eig(Gd-Hd*K) F=Cd*(I-Gd+Hd*K)^-1*Hd Kr=F'*(F*F')^-1
MATLAB nos regresa: K= 61.9203 21.3595
4.1105 -30.2407 -15.6371 -4.8131
-35.8564 -10.6509 -2.8546 41.1216 17.9169 19.7899
6.0904
4.8849
1.5170 -19.7836 -8.7956 -2.4373
ans = 0.8805 + 0.1062i 0.8805 - 0.1062i 0.8024 + 0.1627i 0.8024 - 0.1627i 0.6188 0.6703 F= 0.0273
0.0170 -0.0053
0.0242
0.0352 -0.0133
Kr = 62.7100 -29.5970 -35.5946 41.2589 19.7341 -19.7965
c)
Con el diseño de la parte b) graficaremos el vector de Salida y(kT), el vector de estado x(kT) y el vector de control u(kT) ante un vector de referencia r(kT) de un escalon unitario y un vector de estado de tiempo inicial
El código en Matlab para obtener la Salidas y(kT) Gr=Gd-Hd*K; Hr=Hd*Kr; Cr=Cd; Dr=[0 0;0 0]; X0=[-0.2 0.05 0 0.2 -0.05 0]'
.}
N=100; k=0:N-1; [Y1,X1]=dstep(Gr,Hr,Cr,Dr,1,N); [Y2,X2]=dstep(Gr,Hr,Cr,Dr,2,N); [Y3,X3]=dinitial(Gr,Hr,Cr,Dr,X0,N); Y=(Y1+Y2+Y3); X=(X1+X2+X3); hold off stairs(k*Ts,Y(:,1)) grid on hold on stairs(k*Ts,Y(:,2)) hold on
Matlab nos devuelve el siguiente gráfico:
Para observar los vectores de estado x(kT): Gr=Gd-Hd*K; Hr=Hd*Kr; Cr=Cd; Dr=[0 0;0 0]; X0=[-0.2 0.05 0 0.2 -0.05 0]' N=100;
k=0:N-1; [Y1,X1]=dstep(Gr,Hr,Cr,Dr,1,N); [Y2,X2]=dstep(Gr,Hr,Cr,Dr,2,N); [Y3,X3]=dinitial(Gr,Hr,Cr,Dr,X0,N); Y=(Y1+Y2+Y3); X=(X1+X2+X3); hold off stairs(k*Ts,X(:,1)) grid on hold on stairs(k*Ts,X(:,2)) hold on stairs(k*Ts,X(:,3)) hold on stairs(k*Ts,X(:,4)) hold on stairs(k*Ts,X(:,5)) hold on stairs(k*Ts,X(:,6)) hold on
Matlab nos ofrece el siguiente gráfico:
Para la señal de control U:
La respuesta del sistema usando el K usando la función “place”:
El código en Matlab para obtener la Salidas y(kT) Gr=Gd-Hd*K; Hr=Hd*Kr; Cr=Cd; Dr=[0 0;0 0]; X0=[-0.2 0.05 0 0.2 -0.05 0]' N=100; k=0:N-1; [Y1,X1]=dstep(Gr,Hr,Cr,Dr,1,N); [Y2,X2]=dstep(Gr,Hr,Cr,Dr,2,N); [Y3,X3]=dinitial(Gr,Hr,Cr,Dr,X0,N); Y=(Y1+Y2+Y3); X=(X1+X2+X3); hold off stairs(k*Ts,Y(:,1)) grid on hold on stairs(k*Ts,Y(:,2)) hold on
Matlab nos devuelve el siguiente gráfico:
Para observar los vectores de estado x(kT) : Gr=Gd-Hd*K; Hr=Hd*Kr; Cr=Cd; Dr=[0 0;0 0]; X0=[-0.2 0.05 0 0.2 -0.05 0]' N=100; k=0:N-1; [Y1,X1]=dstep(Gr,Hr,Cr,Dr,1,N); [Y2,X2]=dstep(Gr,Hr,Cr,Dr,2,N); [Y3,X3]=dinitial(Gr,Hr,Cr,Dr,X0,N); Y=(Y1+Y2+Y3); X=(X1+X2+X3); hold off stairs(k*Ts,X(:,1)) grid on hold on stairs(k*Ts,X(:,2)) hold on stairs(k*Ts,X(:,3)) hold on stairs(k*Ts,X(:,4)) hold on
stairs(k*Ts,X(:,5)) hold on stairs(k*Ts,X(:,6)) hold on
Matlab nos ofrece el siguiente gráfico:
Para la señal de control U:
d) Diseñaremos un observador de estados de ganancia Ke usando el método de reubicación de polos para que la dinámica o respuesta de estimación sea lo suficientemente rápida respecto del controlador diseñado en (b). Hallamos en primer lugar la matriz de observabilidad N, lo cual ayudándonos del MATLAB nos da una matriz de orden 6, por tanto la observabilidad del sistema es completa. Sean los polos ubicados en el plano S para el observador:
El código en MATLAB es: %Matriz de Observabilidad N=[Cd;Cd*Gd;Cd*Gd^2;Cd*Gd^3;Cd*Gd^4;Cd*Gd^5] rank(N) se=[-20+4i -20-4i -16+6i -16-6i -30 -40] ze=exp(se*Ts)
MATLAB nos devuelve: >> observador N= 1.0000 0
0 0
0 0
0
0
0
1.0000
0
0
0.9999 0.0400 0.0008 -0.0001 0.0000 -0.0001 -0.0001 0.0000 -0.0001 0.9999 0.0400 0.0008 0.9994 0.0797 0.0029 -0.0004 0.0000 -0.0004 -0.0006 0.0000 -0.0005 0.9996 0.0798 0.0036 0.9982 0.1192 0.0064 -0.0014 0.0000 -0.0015 -0.0019 0.0001 -0.0018 0.9987 0.1193 0.0088 0.9961 0.1580 0.0111 -0.0029 0.0002 -0.0036 -0.0045 0.0003 -0.0043 0.9968 0.1582 0.0169 0.9931 0.1961 0.0172 -0.0053 0.0004 -0.0072 -0.0088 0.0008 -0.0087 0.9939 0.1962 0.0289 ans = 6 se = -20.0000 + 4.0000i -20.0000 - 4.0000i -16.0000 + 6.0000i -16.0000 - 6.0000i -30.0000
-40.0000
ze = 0.4436 + 0.0716i 0.4436 - 0.0716i 0.5122 + 0.1253i 0.5122 - 0.1253i 0.3012
0.2019
Los polos en el plano Z resultan ser:
La ganancia Ke del observador la determinaremos usando el algoritmo del diseño del controlador para sistemas multivariables teniendo en cuenta el criterio de “Sistema Dual”. Para esto es necesario usar las matrices
del sistema original.
La matriz Kes arbitraria para el sistema dual de dimensión 2x6 es:
La matriz
:
Usando el Matlab: Kes=[2 50 5 0 1 -1;1 3 40 0 -1 4] Ge0=Gd'-Cd'*Kes m=eig(G0)
El Matlab no regresa: Kes = 2
50
5
0
1
-1
1
3
40
0
-1
4
Ge0 = -1.0001 -50.0057 -5.2698 -0.0001 -1.0054
0.7328
0.0400 0.9976 -0.1174 0.0000 0.0001 0.0091 0.0008
0.0372
0.8728 -0.0001 -0.0048 -0.2427
-1.0001 -3.0043 -40.2033
0.9999
0.9961 -4.1895
0.0000 0.0000 0.0037 0.0400 0.9983 -0.0912 -0.0001 -0.0040 -0.2023 0.0008 0.0437 1.1971 m= 0.3798 1.2977 + 0.2744i 1.2977 - 0.2744i 0.9625 0.9925 1.0713
Por tanto la Matriz :
[
Los valores propios de la matriz
son distintos, por lo tanto la matriz Kes es v álida. El vector v de
transformación a un sistema de una entrada puede ser:
La matriz
La matriz
numericamente es:
v=[-1 2]' C0=Cd'*v Ke0=place(Ge0,C0,ze) Ke=Ks+v*Ke0
El MATLAB nos regresa: v= -1 2 C0 = -1 0 0 2
de realimentación para el sistema univariable
El código en MATLAB será:
]
es:
0 0 Ke0 = -0.4777
3.3170 -51.6385
0.5867
8.8016 74.2686
Ke = 2.4777 46.6830 56.6385 -0.5867 -7.8016 -75.2686 0.0446
9.6339 -63.2770
1.1734 16.6031 152.5372
Por lo tanto la matriz de ganancia multivariable se obtiene como:
[ ]
del estimador u observador de estados
Aplicando la función “place” del MATLAB se obtiene Ke=place(Gd',Cd',ze)
El MATLAB regresa: Ke = 1.5758 17.3814 51.7124
0.0199 -0.4734 -34.0778
-0.2203 -7.7137 -132.5347
2.0753 33.8543 319.8262
̃
e) Con diseño del observador de estados de la parte d) graficaremos el vector de error entre el vector de estado real
y el vector de estado estimado
un vector de estado real de tiempo inicial el vector de estado estimado inicial nulo.
El Código en MATLAB para los gráficos sería: X0=[-0.2 0.05 0 0.2 -0.05 0]'; N=20; k=0:N-1; e0=X0;
considerando
y
e(:,1)=[0 0 0 0 0 0]; e(:,2)=(Gd-Ke'*Cd)*e0; for i=2:19 e(:,i+1)=(Gd-Ke'*Cd)*e(:,i); end hold off plot(k*Ts,e(1,:)) grid on hold on plot(k*Ts,e(2,:)) hold on plot(k*Ts,e(3,:)) hold on plot(k*Ts,e(4,:)) hold on plot(k*Ts,e(5,:)) hold on plot(k*Ts,e(6,:)) hold on
El Matlab nos regresa el siguiente Gráfico:
You're Reading a Preview Unlock full access with a free trial.
Download With Free Trial
You're Reading a Preview Unlock full access with a free trial.
Download With Free Trial El gráfico que nos da el MATLAB usando la función “place”:
You're Reading a Preview Unlock full access with a free trial.
Download With Free Trial
