CONTROL III EE625-M

La figura 01 describe un sistema mecánico que corresponde al giro de dos bloques inerciales con el
objetivo de analizar las vibraciones en sistemas giratorios y poder controlarlas. Los bloques son
cilindros sólidos de masas M1 , M2, longitudes L1 , L2 y radios R1 y R2 respectivamente.
Los dos bloques son sometidos a torsiones exteriores τ1 y τ2 y desarrollan desplazamientos angulares
θ1 y θ2 con sus respectivas velocidades ω1 y ω2. Ambos bloques están acoplados mediante un resorte
KA y a la vez los bloques están conectados a la base con los resortes K1 y K2 respectivamente.
τ1(t) τ2(t)
B1 B2
K1 KA K2

1(t) R1 R2 2(t)
M1 M2
1(t) 2(t)
L1 L2

Fig 01 Sistema de bloques inerciales.

Las magnitudes de los componentes del sistema están en el sistema MKS y son:
M1 = 4.5 Kg , M2 = 4 Kg
L1 = 0.75 mt , L2 = 0.70 mt
R1 = 0.25 mt , R2 = 0.20mt
B1 = 0.8 N.m.s/rad, B2 = 0.75 N.m.s/rad
K1 =2.5 N.m/rad, K2 = 2.5 N.m/rad , KA =3 N.m/rad

τ1(t) = torsión externa aplicada al bloque 1.

τ2(t) = torsión externa aplicada al bloque 2.
1(t) = desplazamiento angular desarrollado en el eje de la masa 1.
2(t) = desplazamiento angular desarrollado en el eje de la masa 2.
1(t) = velocidad angular desarrollada en el en el eje de la masa 1.
2(t) = velocidad angular desarrollada en el en el eje de la masa 2.

Integrador Método de discretización Número del listado (No)

de la planta
En adelanto zoh Inpar

No = 23
 r(k) e(k) v(k) u(k) Sistema Físico y(k)
Integrador KI

x(k)
K
Km

Fig. 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
KIm = ganancia de realimentación (medición)

r(k) e(k) v(k) u(k) Sistema Físico y(k)

Integrador KI

Estimador
de estados

~
x (k )
K

Km

Fig. 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
KIm = ganancia de realimentación (medición)
 CUARTA PRACTICA CALIFICADA DOMICILIARIA

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. Utilizar una frecuencia de muestreo apropiada para el sistema.

SOLUCION:

La ecuación de estado presenta la siguiente formulación:

𝑥̇ (𝑡) = 𝐴𝑥(𝑡) + 𝐵𝑢(𝑡) … (𝟏)

𝑦(𝑡) = 𝐶𝑥(𝑡) + 𝐷𝑢(𝑡) … (𝟐)

De las ecuaciones (1) y (2) se despeja los valores de 𝜔̇ y 𝜔̇ .

𝜏 (𝐾 + 𝐾 ) 𝐾 𝐵
𝜔̇ = − 𝜃 + 𝜃 − 𝜔 … (𝟑)
𝐽 𝐽 𝐽 𝐽

𝜏 𝐾 (𝐾 + 𝐾 ) 𝐵
𝜔̇ = + 𝜃 − 𝜃 − 𝜔 … (𝟒)
𝐽 𝐽 𝐽 𝐽

Donde:

𝐽 = + = 0.2813
𝑀𝑅 𝑀𝐿
𝐽 = + = 0.2033
4 12

Se realiza cambio de variables:

𝑋 = (𝑥 𝑥 𝑥 𝑥 ) = (𝜃 𝜃 𝜔 𝜔 )

𝑈 = (𝑢 𝑢 ) = (𝜏 𝜏 ) , 𝑦 = (𝑦 𝑦 ) = (𝜃 𝜃 )

Además
𝜔 = 𝜃̇ … (𝟓)

𝜔 = 𝜃̇ … (𝟔)
Expresando en la forma de espacios de estado, se obtiene:

0 0 1 0 0 0
̇𝜃 0 0 0 1 0 0
⎛ (𝐾 + 𝐾 ) 𝐾 𝐵 ⎞ 𝜃 ⎛1 ⎞
𝜏
⎛ 𝜃̇ ⎞ = ⎜− − 0 ⎟ 𝜃 + ⎜𝐽 0⎟
⎜ ⎟ 𝜔
𝜔̇ ⎜ 𝐽 𝐽 𝐽 ⎟ ⎜ ⎟ 𝜏
⎝𝜔̇ ⎠ 𝐾 (𝐾 + 𝐾 ) 𝐵 𝜔 1
− 0 − 0
⎝ 𝐽 𝐽 𝐽 ⎠ ⎝ 𝐽⎠

𝜃
𝜃 1 0 0 0 𝜃
=
𝜃 0 1 0 0 𝜔
𝜔
𝜃̇ 0 0 1 0 𝜃
⎛ 𝜃̇ ⎞ = 0 0 0 1 𝜃
𝜔̇ −19.5521 10.6648 − 2.8439 0 𝜔
⎝𝜔̇ ⎠ 14.7565 −27.0536 0 −3.6891 𝜔

0 0
0 0 𝜏
+ 𝜏 … (𝟕)
3.5549 0
0 4.9188
𝜃
𝜃 1 0 0 0 𝜃
= … (𝟖)
𝜃 0 1 0 0 𝜔
𝜔

Luego se emplea el cambio de variable:

𝑥̇ (𝑡) 0 0 1 0 𝑥 (𝑡)
⎛𝑥̇ (𝑡)⎞ 0 0 0 1 ⎛𝑥 (𝑡)⎞
=
𝑥̇ (𝑡) −19.5521 10.6648 − 2.8439 0 𝑥 (𝑡)
⎝𝑥̇ (𝑡)⎠ 14.7565 −27.0536 0 −3.6891 ⎝𝑥 (𝑡)⎠

0 0
0 0 𝑢
+ 𝑢
3.5549 0
0 4.9188

𝑥 (𝑡)
𝑦 1 0 0 0 ⎛𝑥 (𝑡)⎞ 0 0 𝑢
𝑦 = + 𝑢
0 1 0 0 𝑥 (𝑡) 0 0
⎝𝑥 (𝑡)⎠
Diagrama de bloques realizado en Simulink:

b) Determinar la Matriz de Transferencia Pulso del sistema en lazo abierto.

SOLUCION:

De la Ecuación de Estados:
𝑥̇ = 𝐴𝑥 + 𝐵𝑢
𝑦 = 𝐶𝑥 + 𝐷𝑢
Pasando al dominio de Laplace:

𝑠𝑥 = 𝐴𝑥 + 𝐵𝑢 → (𝑠𝐼 − 𝐴)𝑥 = 𝐵𝑢 → 𝑥 = (𝑠𝐼 − 𝐴) 𝐵𝑢

𝑦 = 𝐶𝑥 + 𝐷𝑢 → 𝑦 = (𝐶(𝑠𝐼 − 𝐴) 𝐵 + 𝐷)𝑢

Siendo la matriz de transferencia:

(𝐶(𝑠𝐼 − 𝐴) 𝐵 + 𝐷)

Se emplea el siguiente código en MATLAB:

clc
clear
close all
syms s
A=[ 0 0 1 0
0 0 0 1
-19.5521 10.6648 -2.8439 0
14.7565 -27.0536 0 -3.6891];

B=[0 0
0 0
3.5549 0
0 4.9188];
C=[1 0 0 0
0 1 0 0];
D=[0 0; 0 0];
[Nuz, Dez]=ss2tf(A,B,C,D,1)
[Nuz, Dez]=ss2tf(A,B,C,D,2)

Se obtiene:

3.5549𝑠 + 13.1144𝑠 + 96.1728 52.4580

𝜃 52.4579 4.9188𝑠 + 13.9886𝑠 + 96.1729 𝜏
= 𝜏
𝜃 𝑠 + 6.5330𝑠 + 57.0971𝑠 + 149.0674𝑠 + 371.5796

c) Diseñar un sistema de control de seguimiento con realimentación de estado de ganancia K y

un integrador de ganancia KI usando Control Optimo Cuadrático en tiempo discreto según
la figura (02), eligiendo adecuadamente las matrices de ponderación Q y R y la matriz Km
considerando que el rango del vector r(kT) varía desde -12V a +12V, mientras que el rango
de las variables de salida y(kT) varía desde -180° a +180° en estado estacionario,

SOLUCION:

Las matrices ampliadas debido al efecto de integrador ya fueron calculadas y se muestran:

1 0 0 0 0 0
𝐶 =
0 1 0 0 0 0

0 0
𝐷 =( )
0 0
Se comprueba la controlabilidad y rango con el código en MATLAB:

Ma=ctrb(Ga,Ha
size(Ma)
rM=rank(Ma
La solución de la ecuación de RICATTI se obtiene mediante la ecuación:
𝑃(𝑘 + 1) = 𝑄 + 𝐺 𝑃(𝑘)𝐺 − 𝐺 𝑃(𝑘)𝐻 𝑅 + 𝐻 𝑃(𝑘 )𝐻 𝐻 𝑃(𝑘)𝐺 … (1)

Elegimos las matrices de ponderación 𝑄 𝑦 𝑅:

100 0 0 0 0 0
0 100 0 0 0 0
⎛ 0 0 90 0 0 0 ⎞ 0.25 0
𝑄 =⎜ ⎟, 𝑅 =
⎜ 0 0 0 90 0 0 ⎟ 0 0.25
0 0 0 0 150 0
⎝ 0 0 0 0 0 150 ⎠
Entonces procedemos a los cálculos numéricos, para esto nos basamos en las ecuaciones de estado y
salida en tiempo discreto de la planta obtenidas y aplicamos algorítmicamente la ecuación (1),
mediante el siguiente código elaborado en el programa MATLAB(se adjunta en ANEXOS):

%% Diseño de un controlador optimo con integrador

km=3.8197;
Ga=[G [0;0;0;0] [0;0;0;0];-km*C*G [1;0] [0;1]]
Ha=[H;-km*C*H]
Ca=[C [0;0] [0;0]]
Da=D
Qa=diag([100 100 90 90 150 150]);
R = [0.25 0;0 0.25];
P = zeros(6);
n=10;
for i = 1:n % Iteraciones se revisa cada 10 o 20 iteraciones
P = Qa + Ga'*P*Ga - Ga'*P*Ha*inv(R + Ha'*P*Ha)*Ha'*P*Ga
end;
% Matriz de Ganancia de Control Óptimo
K1 = inv(R + Ha'*P*Ha)*Ha'*P*Ga
K = K1(:,1:4)
KI = -K1(:,5:6)
Se obtiene:

A través de MATLAB se obtiene la respuesta del sistema ante una entrada de escalón:

r=1;
[Yr Xar]=dstep(Gac,Hac*r,Cac,Dac,1,Nm);
k=0:Nm-1;
figure
stairs(k*Ts,Yr)
title({'Respuesta del sistema ante entrada de escalon','CONTROL OPTIMO'})
ylabel('Yr(kT)')
xlabel('kT(seg)')
grid
Variando la matriz Qa de acuerdo la dinámica del sistema se presenta los siguientes resultados:

Un integrador más lento

100 0 0 0 0 0
0 100 0 0 0 0
⎛ 0 0 90 0 0 0 ⎞ 0.25 0
𝑄 =⎜ ⎟, 𝑅 =
⎜ 0 0 0 90 0 0 ⎟ 0 0.25
0 0 0 0 50 0
⎝ 0 0 0 0 0 50 ⎠
Se observa un mayor sobreimpulso

Un integrador más fino

100 0 0 0 0 0
0 100 0 0 0 0
⎛ 0 0 90 0 0 0 ⎞ 0.25 0
𝑄 = 0.01 ⎜ ⎟, 𝑅 =
⎜ 0 0 0 90 0 0 ⎟ 0 0.25
0 0 0 0 50 0
⎝ 0 0 0 0 0 50 ⎠
Se observa un menor sobreimpulso y más lento en estabilizarse.

Un integrador más rápido

100 0 0 0 0 0
0 100 0 0 0 0
⎛ 0 0 90 0 0 0 ⎞ 0.25 0
𝑄 = 0.01 ⎜ ⎟, 𝑅 =
⎜ 0 0 0 90 0 0 ⎟ 0 0.25
0 0 0 0 1500 0
⎝ 0 0 0 0 0 1500 ⎠
Se observa muy menor sobreimpulso.

d) Determinar las funciones de transferencia de lazo cerrado,

SOLUCION:

La matriz de transferencia de lazo cerrado es:

𝑌(𝑧)
𝐹(𝑧) = = 𝐶 (𝑧𝐼 − 𝐺 ) 𝐻 +𝐷
𝑈(𝑧)
Donde:

Se emplea el siguiente código de MATLAB para determinar la función de transferencia de lazo

cerrado:

%%Funcion de transferencia lazo cerrado

Gac=Ga-Ha*K1
Hac=Ha
Cac=Ca-Da*K1
Dac=zeros(2,2)
M=ctrb(Gac,Hac);%verificación de la controlabilidad
rM=rank(M)
[Nuz, Dez]=ss2tf(Gac,Hac,Cac,Dac,1)
[Nuz, Dez]=ss2tf(Gac,Hac,Cac,Dac,2)
Se reemplaza las matrices en la expresión de 𝑭(𝒛):

𝑀 𝑃
𝑁 𝑄
𝑧 −0.0887𝑧 − 0.03192𝑧 + 0.006344𝑧 − 0.000246𝑧 − 0.0000021𝑧 + 0.0000000065

Donde:

𝑀 = 0.29678𝑧 − 0.29986𝑧 − 0.001608𝑧 + 0.01211𝑧 − 0.00133z + 0.000019

𝑃 = 0.17809𝑧 − 0.18708𝑧 + 0.00245𝑧 + 0.00773𝑧 − 0.001207z + 0.000018
𝑁 = 0.17809𝑧 − 0.18549𝑧 + 0.0002381𝑧 + 0.008427𝑧 − 0.00128z + 0.0000016
𝑄 = 0.26977𝑧 − 0.27500𝑧 − 0.006175𝑧 + 0.01282𝑧 − 0.001432z + 0.0000015

e) Graficar las variables de estado x(kT), del error e(kT), v(k) y(kT) y el vector de control u(kT)
ante un vector de entrada en r(kT) = [8 10]TV (tipo escalón) usando como vector de
condiciones iniciales x(0)  0 según su buen criterio

SOLUCION:
Figuras (15) a la (19)

Se considera la siguiente condición inicial: x0=[30 30 6 6 0 0]'

El código realizado en MATLAB se adjunta en los ANEXOS

Graficas Obtenidas:

VECTOR DE ESTADOS
 VARIABLE DE ERROR

VARIABLE ACUMULADOR DE ERROR

 VECTOR DE SALIDA

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

Al variar la matriz variable de ponderación Qa de acuerdo a la dinámica del sistema se puede

interpretar físicamente de siguiente manera:
Para un integrador más lento se obtiene un sobreimpulso mayor.
Para un integrador más fino se obtiene un menor sobreimpulso y más lento para estabilizarse.
Para un integrador más rápido se obtiene muy menor sobreimpulso.
Se puede apreciar que aproximadamente 7 segundos el sistema sigue prácticamente a la referencia
y error tiende a cero.
La matriz de prealimentación hace que las salidas sigan a la entrada y a la vez transforma las
entradas de referencia en entradas de control.
Las variables de estado alcanzan valores elevados en la zona transitoria los cuales se pueden
reducir reubicando los polos.

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, de tal manera que el
error e(k )  x(k )  ~
x (k ) tienda a cero a una velocidad igual a 3 veces o mayor que la velocidad
de respuesta de la salida del sistema, mediante la técnica de ubicación de polos para el diseño
del estimador Ke usando el método algorítmico

SOLUCION:

La Ganancia Ke del observador se determina utilizando 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 𝐺 y 𝐶 del sistema original La matriz Kes arbitraria para el sistema dual de dimensión
2×4 es:

Los valores propios de la matriz Ge0 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 Ke0 de realimentación para el sistema univariable (Ge0, He0) es:

Se utilizo el siguiente código en el programa MATLAB:

%% ESTIMADOR O OBSERVADOR de estados (g)

Kes=[-1 2 0 -2
-5 2 0 4 ];
ve=[-1; -2];
Ge0 = G'-C'*Kes
He0 = C'*ve
%el Ke0 para el sistema univariale Ge0, He0
M=[He0 Ge0*He0 Ge0^2*He0 Ge0^3*He0];
S = [-0.801+1.1847i, -0.801-1.1847i, -5, -5];%velocidad 3 veces o mayor que el sistema
rM=rank(M) ;
Pa=poly(Ge0);
a1=Pa(2); a2=Pa(3); a3=Pa(4); a4=Pa(5);
W=[a3 a2 a1 1;a2 a1 1 0;a1 1 0 0;1 0 0 0];
T=M*W;
Z=exp(Ts*S);
Pc=poly(Z);
al1=Pc(2); al2=Pc(3); al3=Pc(4); al4=Pc(5);
Ke0=[al4-a4 al3-a3 al2-a2 al1-a1]*inv(T)
Ke=(Kes+ve*Ke0)'
h) Graficar las variables de estado x(k), ~x (t ) , del error e(kT), v(k) y(kT) y el vector de control
u(kT) ante un vector de entrada en r(kT) = [8 10]TV (tipo escalón) usando como vector de
condiciones iniciales e(0)=(-6 8 5 -4)T.
SOLUCION:

Figuras (20) a la (25)

Se considera las siguientes condiciones iniciales:

x0=[30 30 6 6]';

e0=[-6 8 5 -4]';

El código realizado en MATLAB se adjunta en los ANEXOS

Graficas Obtenidas:
r(kT) = 8 V (tipo escalón)

VECTOR DE ESTADOS
VECTOR ESTIMADOR DE ESTADOS

VARIABLE DE ERROR
VARIABLE ACUMULADOR DE ERROR

VECTOR DE SALIDA
 VARIABLE DE CONTROL

r(kT) = 10 V (tipo escalón)

VECTOR DE ESTADOS
VECTOR ESTIMADOR DE ESTADOS

VARIABLE DE ERROR
VARIABLE ACUMULADOR DE ERROR

VECTOR DE SALIDA
 VARIABLE DE CONTROL

i) Comparar e interpretar físicamente la respuesta del sistema basándose en las gráficas de las
partes (e) y (h).

El diseño de la ubicación de polos para la planta es independiente del diseño del observador, lo
cual permite dividir el proceso de diseño en dos etapas separadas, primero se determina la matriz
de ganancia de realimentación K para obtener la ubicación de polos deseada para la planta, y a
continuación se determinan la matriz Ko para la ubicación de polos deseados del observador.
Se deduce que el número total de polos del conjunto planta – observador de estado es la suma de
los productos por cada uno por separado, con lo cual, si los correspondiente al observador están
situados lo suficientemente a la izquierda de los correspondientes a la planta, su respuesta
transitoria decaerá mucho más rápidamente y, como consecuencia, los polos dominantes serán los
correspondientes a la planta.
La ventaja de la estimación de estado es que evita la necesidad de utilizar varios sensores para
medir las variables de estado, si estas fueran accesibles.

j) Determinar la Función de Transferencia Pulso del controlador-estimador U(z)/-Y(z)

SOLUCION:

La matriz de transferencia de lazo cerrado es:

𝑌(𝑧)
𝐹(𝑧) = = 𝐶 (𝑧𝐼 − 𝐺 ) 𝐻 +𝐷
𝑈(𝑧)
Donde:

Se emplea el siguiente código de MATLAB para determinar la función de transferencia de lazo

cerrado:

%%Funcion de transferencia pulso del controlador-estimador

M=ctrb(Ge0,He0);%verificación de la controlabilidad
rM=rank(M)
[Nuz, Dez]=ss2tf(Ge0,He0,Cec,Dec,1)

Se reemplaza las matrices en la expresión de 𝑭(𝒛):

𝑀
𝑁
𝑧 + 1.0614𝑧 + 8.8002𝑧 + 0.0412𝑧 − 0.3466

Donde:

𝑀 = −𝑧 + 2.3435𝑧 + 0.0985z − 0.0771

𝑁 = −2𝑧 − 3.1163𝑧 + 0.0795z + 0.1361
 ANEXOS

CODIGO GENERAL COMPLETO 3PRACTICA Y 4PRACTICA

clc
clear
close all
%% Datos
% (a)
A=[ 0 0 1 0
0 0 0 1
-19.5521 10.6648 -2.8439 0
14.7565 -27.0536 0 -3.6891];
B=[0 0
0 0
3.5549 0
0 4.9188];
C=[1 0 0 0
0 1 0 0];
D=[0 0; 0 0];
%t=0:0.1:15;
%step(A,B,C,D,2,15)
%grid
Ts=1;
[G H C D]=c2dm(A,B,C,D,Ts,'zoh')
%M=ctrb(G,H) %verificación de la controlabilidad
%rM=rank(M)
%[Nuz, Dez]=ss2tf(G,H,C,D,1)
Ks=[-1 4 0 -2
-5 4 0 4 ];
v=[-1; -4];
G0=G-H*Ks
H0=H*v
M=[H0 G0*H0 G0^2*H0 G0^3*H0]
S = [-0.801+1.1847i, -0.801-1.1847i, -3, -4];
rM=rank(M)
Pa=poly(G0)
a1=Pa(2); a2=Pa(3); a3=Pa(4); a4=Pa(5);
W=[a3 a2 a1 1;a2 a1 1 0;a1 1 0 0;1 0 0 0]
T=M*W
Z=exp(Ts*S)
Pc=poly(Z)
al1=Pc(2); al2=Pc(3); al3=Pc(4); al4=Pc(5);
K0=[al4-a4 al3-a3 al2-a2 al1-a1]*inv(T)
K=Ks+v*K0
F=C*inv(eye(4)-G+H*K)*H
Kr=inv(F)
%% Simulacion como si fuera un sistema SISO
%diseño de regulador con realimentacion de estado de ganancia K
Gc=G-H*K%matriz de lazo cerrado
Hc=H*Kr
Cc=C-D*K
Dc=D*Kr
figure
dstep(Gc,Hc,Cc,Dc)% respuesta ante un escalon
Nm=10;
[Yd1,Xd1]=dstep(Gc,Hc,Cc,Dc,1,Nm);
[Yd2,Xd2]=dstep(Gc,Hc,Cc,Dc,2,Nm);
Y=Yd1+Yd2;
X=Xd1+Xd2;
k=0:Nm-1;
figure
stairs(k*Ts,Y)
figure
stairs(k*Ts,X)
%agregando integradores con la funcion place para sistema MIMO
%Sistema de seguimiento incluyendo integradores
km=3.8197;
Ga=[G [0;0;0;0] [0;0;0;0];-km*C*G [1;0] [0;1]]
Ha=[H;-km*C*H]
Ca=[C [0;0] [0;0]]
Da=D
Ma=ctrb(Ga,Ha)
size(Ma)
rM=rank(Ma)
Sa = [-0.801+1.1847i, -0.801-1.1847i, -3, -4, -5, -6];
Za=exp(Ts*Sa);
Ka=place(Ga,Ha,Za)
K=Ka(:,1:4)
Ki=-Ka(:,5:6)
Gac=Ga-Ha*Ka
Hac=Ha
Cac=Ca-Da*Ka
Dac=zeros(2,2)
M=ctrb(Gac,Hac) %verificación de la controlabilidad
rM=rank(M)
%[Nuz, Dez]=ss2tf(Gac,Hac,Cac,Dac,2)
dstep(Gac,Hac,Cac,Dac)
Nm=10;
%(h) con condiciones iniciales
x0=[30 30 6 6 0 0]';
r1=8;
r2=10;
r1=(r1*ones(size(k)));
r2=(r2*ones(size(k)));
R=[r1;r2];
[Y1,X1]=dlsim(Gac,Hac,Cac,Dac,R,x0);%ante cualquier entrada arbitraria incluyendo ci
[Y2,X2]=dlsim(Gac,Hac,Cac,Dac,R,x0);
Km=3.8197;
Y=Y1+Y2;
X=X1+X2;
E=R'-Km*Y;
U=(-Ka*X')';
V=[X1(:,6),X2(:,6)];%V Acumulacion del error ante entrada de escalon
figure
stairs(k*Ts,X)
title({'Respuesta del sistema ante entrada r(KT)=[8;10]V','Vector de estados'})
ylabel('Vector de estados X(kT)')
xlabel('kT(seg)')
legend('x_1','x_2','x_3','x_4')
grid
figure
stairs(k*Ts,E)
title({'Respuesta del sistema ante entrada r(KT)=[8;10] V','Variable de error'})
ylabel('Variable de error E(kT)')
xlabel('kT(seg)')
grid
figure
stairs(k*Ts,V)
title({'Respuesta del sistema ante entrada r(KT)=[8;10] V','Variable acumulacion de error'})
ylabel('Variable de acumulacion del error V(kT)')
xlabel('kT(seg)')
grid
figure
stairs(k*Ts,Y)
title({'Respuesta del sistema ante entrada r(KT)=[8;10] V','Vector de Salida'})
ylabel('Y(kT)')
xlabel('kT(seg)')
grid
figure
stairs(k*Ts,U)
title({'Respuesta del sistema ante entrada r(KT)=[8;10] V','Variable de control'})
ylabel('Variable de control U(kT)')
xlabel('kT(seg)')
grid
%(i) condiciones iniciales nulas
Nm=20;
x0=[0 0 0 0 0 0]';
r1=7.5;
r2=7.5;
k=0:25-1;
R1=r1*ones(size(k));
R2=r2*ones(size(k));
R1=r1.*sin(1.25*Ts*k);%variar frecuencia
R2=r2.*sign(cos(0.79*Ts*k));
R=[R1;R2];
[Y1,X1]=dlsim(Gac,Hac,Cac,Dac,R,x0);
[Y2,X2]=dlsim(Gac,Hac,Cac,Dac,R,x0);
Km=3.8197;
Y=Y1+Y2;
X=X1+X2;
E=R'-Km*Y;
U=(-Ka*X')';
V=[X1(:,6),X2(:,6)];%V Acumulacion del error ante entrada de escalon
figure
stairs(k*Ts,X)
title({'Respuesta del sistema ante entrada r(KT)=[7.5sinwt 7.5sign(coswt)]','Vector de
estados'})
ylabel('Vector de estados X(kT)')
xlabel('kT(seg)')
legend('x_1','x_2','x_3','x_4')
grid
figure
stairs(k*Ts,E)
title({'Respuesta del sistema ante entrada r(KT)=[7.5sinwt 7.5sign(coswt)]','Variable de
error'})
ylabel('Variable de error E(kT)')
xlabel('kT(seg)')
grid
figure
stairs(k*Ts,V)
title({'Respuesta del sistema ante entrada r(KT)=[7.5sinwt 7.5sign(coswt)]','Variable
acumulacion de error'})
ylabel('Variable de acumulacion del error V(kT)')
xlabel('kT(seg)')
grid
figure
stairs(k*Ts,Y)
title({'Respuesta del sistema ante entrada r(KT)=[7.5sinwt 7.5sign(coswt)]','Vector de
Salida'})
ylabel('Y(kT)')
xlabel('kT(seg)')
grid
figure
stairs(k*Ts,U)
title({'Respuesta del sistema ante entrada r(KT)=[7.5sinwt 7.5sign(coswt)]','Variable de
control'})
ylabel('Variable de control U(kT)')
xlabel('kT(seg)')
grid
%% Diseño de un controlador OPTIMO con integrador
km=3.8197;
Ga=[G [0;0;0;0] [0;0;0;0];-km*C*G [1;0] [0;1]]
Ha=[H;-km*C*H]
Ca=[C [0;0] [0;0]]
Da=D
Qa=diag([100 100 90 90 150 150]);
R = [0.25 0;0 0.25];
P = zeros(6);
n=10;
for i = 1:n % Iteraciones se revisa cada 10 o 20 iteraciones
P = Qa + Ga'*P*Ga - Ga'*P*Ha*inv(R + Ha'*P*Ha)*Ha'*P*Ga
end;
% Matriz de Ganancia de Control Óptimo
K1 = inv(R + Ha'*P*Ha)*Ha'*P*Ga
K = K1(:,1:4)
KI = -K1(:,5:6)
%%Funcion de transferencia lazo cerrado
Gac=Ga-Ha*K1
Hac=Ha
Cac=Ca-Da*K1
Dac=zeros(2,2)
M=ctrb(Gac,Hac);%verificación de la controlabilidad
rM=rank(M)
[Nuz, Dez]=ss2tf(Gac,Hac,Cac,Dac,2);
%Grafica del sistema ante entrada de escalon
r=1;
[Yr Xar]=dstep(Gac,Hac*r,Cac,Dac,1,Nm);
k=0:Nm-1;
figure
stairs(k*Ts,Yr)
title({'Respuesta del sistema ante entrada de escalon','CONTROL OPTIMO'})
ylabel('Yr(kT)')
xlabel('kT(seg)')
grid
%(e) con condiciones iniciales
x0=[30 30 6 6 0 0]';
r1=8;
r2=10;
r1=(r1*ones(size(k)));
r2=(r2*ones(size(k)));
R=[r1;r2];
[Y1,X1]=dlsim(Gac,Hac,Cac,Dac,R,x0);%ante cualquier entrada arbitraria incluyendo ci
[Y2,X2]=dlsim(Gac,Hac,Cac,Dac,R,x0);
Km=3.8197;
Y=Y1+Y2;
X=X1+X2;
E=R'-Km*Y;
U=(-Ka*X')';
V=[X1(:,6),X2(:,6)];%V Acumulacion del error ante entrada de escalon
figure
stairs(k*Ts,X)
title({'Respuesta del sistema ante entrada r(KT)=[8;10]V','Vector de estados','CONTROL
OPTIMO'})
ylabel('Vector de estados X(kT)')
xlabel('kT(seg)')
legend('x_1','x_2','x_3','x_4')
grid
figure
stairs(k*Ts,E)
title({'Respuesta del sistema ante entrada r(KT)=[8;10] V','Variable de error','CONTROL
OPTIMO'})
ylabel('Variable de error E(kT)')
xlabel('kT(seg)')
grid
figure
stairs(k*Ts,V)
title({'Respuesta del sistema ante entrada r(KT)=[8;10] V','Variable acumulacion de
error','CONTROL OPTIMO'})
ylabel('Variable de acumulacion del error V(kT)')
xlabel('kT(seg)')
grid
figure
stairs(k*Ts,Y)
title({'Respuesta del sistema ante entrada r(KT)=[8;10] V','Vector de Salida','CONTROL
OPTIMO'})
ylabel('Y(kT)')
xlabel('kT(seg)')
grid
figure
stairs(k*Ts,U)
title({'Respuesta del sistema ante entrada r(KT)=[8;10] V','Variable de control','CONTROL
OPTIMO'})
ylabel('Variable de control U(kT)')
xlabel('kT(seg)')
grid
%% ESTIMADOR O OBSERVADOR de estados (g)
Kes=[-1 2 0 -2
-5 2 0 4 ];
ve=[-1; -2];
Ge0 = G'-C'*Kes
He0 = C'*ve
%el Ke0 para el sistema univariable Ge0, He0
M=[He0 Ge0*He0 Ge0^2*He0 Ge0^3*He0];
S = [-0.801+1.1847i, -0.801-1.1847i, -5, -5];%velocidad 3 veces o mayor que el sistema
rM=rank(M) ;
Pa=poly(Ge0);
a1=Pa(2); a2=Pa(3); a3=Pa(4); a4=Pa(5);
W=[a3 a2 a1 1;a2 a1 1 0;a1 1 0 0;1 0 0 0];
T=M*W;
Z=exp(Ts*S);
Pc=poly(Z);
al1=Pc(2); al2=Pc(3); al3=Pc(4); al4=Pc(5);
Ke0=[al4-a4 al3-a3 al2-a2 al1-a1]*inv(T)
Ke=(Kes+ve*Ke0)'
%% Simulacion del observador sin realimentar de lazo cerrado
Nm=20;
Gec=G-Ke*C;
Hec=[0;0;0;0];%no hay entrada es regulador
Cec=C;
Dec=[0;0];
%con condiciones iniciales
x0=[30 30 6 6]';
xe=[-6 8 5 -4 ]';%e estimador dato del problema
e0=x0-xe;
[Y0 X0]=dinitial(Gec,Hec*r1,Cec,Dec*r1,x0,Nm);%se mutiplica por 0 porque no hay
entrada solo con CI
V=[X0(:,4)];
U=-Ke0*X0';
r1=10;% primera referencia
figure
stairs(k*Ts,X0)
legend('x1(k)','x2(k)','x3(k)','x4(k)')
title({'Vector de estados del sistema r(KT)=10V','CON ESTIMADOR DE ESTADOS'})
ylabel('X(kT)')
xlabel('kT(seg)')
grid
[Ye Xe]=dinitial(Gec,Hec*r1,Cec,Dec*r1,xe,Nm);
figure
stairs(k*Ts,Xe)
legend('~x1(k)','~x2(k)','~x3(k)','~x4(k)')
title({'Vector de estados estimados del sistema r(KT)=10V','CON ESTIMADOR DE
ESTADOS'})
ylabel('X(kT)')
xlabel('kT(seg)')
grid
[~,E]=dinitial(Gec,Hec*r1,Cec,Dec,e0,Nm);
figure
stairs(k*Ts,E)%salida vector error de estado comparar con el integrar tiempo
establecimiento
title({'Vector de error r(KT)=10V','CON ESTIMADOR DE ESTADOS'})
grid
figure
stairs(k*Ts,V)
title({'Variable acumulacion de error r(KT)=10V','CON ESTIMADOR DE ESTADOS'})
ylabel('Variable de acumulacion del error V(kT)')
xlabel('kT(seg)')
grid
figure
stairs(k*Ts,Y0)
title({'Vector de salida r(KT)=10V','CON ESTIMADOR DE ESTADOS'})
ylabel('Y(kT)')
xlabel('kT(seg)')
grid
figure
stairs(k*Ts,U)
title({'Variable de control r(KT)=10V','CON ESTIMADOR DE ESTADOS'})
ylabel('Variable de control U(kT)')
xlabel('kT(seg)')
grid
%%Funcion de transferencia pulso del controlador-estimador
M=ctrb(Ge0,He0);%verificación de la controlabilidad
rM=rank(M)
[Nuz, Dez]=ss2tf(Ge0,He0,Cec,Dec,1)
 CODIGOS POR PREGUNTAS 4 PRACTICA

PREGUNTA B

%cálculo de matriz de transferencia a lazo abierto

clc
clear
close all
syms s
A=[ 0 0 1 0
0 0 0 1
-19.5521 10.6648 -2.8439 0
14.7565 -27.0536 0 -3.6891];
B=[0 0
0 0
3.5549 0
0 4.9188];
C=[1 0 0 0
0 1 0 0];
D=[0 0; 0 0];

[Nuz, Dez]=ss2tf(A,B,C,D,1)
[Nuz, Dez]=ss2tf(A,B,C,D,2)

PREGUNTA C

%% Diseño de un controlador optimo con integrador

km=3.8197;
Ga=[G [0;0;0;0] [0;0;0;0];-km*C*G [1;0] [0;1]]
Ha=[H;-km*C*H]
Ca=[C [0;0] [0;0]]
Da=D
Qa=diag([100 100 90 90 150 150]);
R = [0.25 0;0 0.25];
P = zeros(6);
n=10;
for i = 1:n % Iteraciones se revisa cada 10 o 20 iteraciones
P = Qa + Ga'*P*Ga - Ga'*P*Ha*inv(R + Ha'*P*Ha)*Ha'*P*Ga
end;
% Matriz de Ganancia de Control Óptimo
K1 = inv(R + Ha'*P*Ha)*Ha'*P*Ga
K = K1(:,1:4)
KI = -K1(:,5:6)

PREGUNTA D

%%Funcion de transferencia lazo cerrado

Gac=Ga-Ha*K1
Hac=Ha
Cac=Ca-Da*K1
Dac=zeros(2,2)
M=ctrb(Gac,Hac);%verificación de la controlabilidad
rM=rank(M)
[Nuz, Dez]=ss2tf(Gac,Hac,Cac,Dac,1)
[Nuz, Dez]=ss2tf(Gac,Hac,Cac,Dac,2)

PROBLEMA E

%(e) con condiciones iniciales

x0=[30 30 6 6 0 0]';
r1=8;
r2=10;
r1=(r1*ones(size(k)));
r2=(r2*ones(size(k)));
R=[r1;r2];
[Y1,X1]=dlsim(Gac,Hac,Cac,Dac,R,x0);%ante cualquier entrada arbitraria incluyendo ci
[Y2,X2]=dlsim(Gac,Hac,Cac,Dac,R,x0);
Km=3.8197;
Y=Y1+Y2;
X=X1+X2;
E=R'-Km*Y;
U=(-Ka*X')';
V=[X1(:,6),X2(:,6)];%V Acumulacion del error ante entrada de escalon
figure
stairs(k*Ts,X)
title({'Respuesta del sistema ante entrada r(KT)=[8;10]V','Vector de estados','CONTROL
OPTIMO'})
ylabel('Vector de estados X(kT)')
xlabel('kT(seg)')
legend('x_1','x_2','x_3','x_4')
grid
figure
stairs(k*Ts,E)
title({'Respuesta del sistema ante entrada r(KT)=[8;10] V','Variable de error','CONTROL
OPTIMO'})
ylabel('Variable de error E(kT)')
xlabel('kT(seg)')
grid
figure
stairs(k*Ts,V)
title({'Respuesta del sistema ante entrada r(KT)=[8;10] V','Variable acumulacion de
error','CONTROL OPTIMO'})
ylabel('Variable de acumulacion del error V(kT)')
xlabel('kT(seg)')
grid
figure
stairs(k*Ts,Y)
title({'Respuesta del sistema ante entrada r(KT)=[8;10] V','Vector de Salida','CONTROL
OPTIMO'})
ylabel('Y(kT)')
xlabel('kT(seg)')
grid
figure
stairs(k*Ts,U)
title({'Respuesta del sistema ante entrada r(KT)=[8;10] V','Variable de control','CONTROL
OPTIMO'})
ylabel('Variable de control U(kT)')
xlabel('kT(seg)')
grid

PROBLEMAS G
%% ESTIMADOR O OBSERVADOR de estados (g)
Kes=[-1 2 0 -2
-5 2 0 4 ];
ve=[-1; -2];
Ge0 = G'-C'*Kes
He0 = C'*ve
%el Ke0 para el sistema univariale Ge0, He0
M=[He0 Ge0*He0 Ge0^2*He0 Ge0^3*He0];
S = [-0.801+1.1847i, -0.801-1.1847i, -5, -5];%velocidad 3 veces o mayor que el sistema
rM=rank(M) ;
Pa=poly(Ge0);
a1=Pa(2); a2=Pa(3); a3=Pa(4); a4=Pa(5);
W=[a3 a2 a1 1;a2 a1 1 0;a1 1 0 0;1 0 0 0];
T=M*W;
Z=exp(Ts*S);
Pc=poly(Z);
al1=Pc(2); al2=Pc(3); al3=Pc(4); al4=Pc(5);
Ke0=[al4-a4 al3-a3 al2-a2 al1-a1]*inv(T)
Ke=(Kes+ve*Ke0)'

PREGUNTA H
r(kT) = 8 V (tipo escalón)

%% Simulacion del observador sin realimentar de lazo cerrado

Nm=20;
Gec=G-Ke*C;
Hec=[0;0;0;0];%no hay entrada es regulador
Cec=C;
Dec=[0;0];
%con condiciones iniciales
x0=[30 30 6 6]';
xe=[-6 8 5 -4 ]';%e estimador dato del problema
e0=x0-xe;
[Y0 X0]=dinitial(Gec,Hec*r1,Cec,Dec*r1,x0,Nm);%se mutiplica por 0 porque no hay
entrada solo con CI
V=[X0(:,4)];
U=-Ke0*X0';
r1=8;% primera referencia
figure
stairs(k*Ts,X0)
legend('x1(k)','x2(k)','x3(k)','x4(k)')
title({'Vector de estados del sistema','CON ESTIMADOR DE ESTADOS'})
ylabel('X(kT)')
xlabel('kT(seg)')
grid
[Ye Xe]=dinitial(Gec,Hec*r1,Cec,Dec*r1,xe,Nm);
figure
stairs(k*Ts,Xe)
legend('~x1(k)','~x2(k)','~x3(k)','~x4(k)')
title({'Vector de estados estimados del sistema','CON ESTIMADOR DE ESTADOS'})
ylabel('X(kT)')
xlabel('kT(seg)')
grid
[~,E]=dinitial(Gec,Hec*r1,Cec,Dec,e0,Nm);
figure
stairs(k*Ts,E)%salida vector error de estado comparar con el integrar tiempo
establecimiento
title({'Vector de error','CON ESTIMADOR DE ESTADOS'})
grid
figure
stairs(k*Ts,V)
title({'Variable acumulacion de error','CON ESTIMADOR DE ESTADOS'})
ylabel('Variable de acumulacion del error V(kT)')
xlabel('kT(seg)')
grid
figure
stairs(k*Ts,Y0)
title({'Vector de salida','CON ESTIMADOR DE ESTADOS'})
ylabel('Y(kT)')
xlabel('kT(seg)')
grid
figure
stairs(k*Ts,U)
title({'Variable de control','CON ESTIMADOR DE ESTADOS'})
ylabel('Variable de control U(kT)')
xlabel('kT(seg)')
grid

r(kT) = 10 V (tipo escalón)

%% Simulacion del observador sin realimentar de lazo cerrado

Nm=20;
Gec=G-Ke*C;
Hec=[0;0;0;0];%no hay entrada es regulador
Cec=C;
Dec=[0;0];
%con condiciones iniciales
x0=[30 30 6 6]';
xe=[-6 8 5 -4 ]';%e estimador dato del problema
e0=x0-xe;
[Y0 X0]=dinitial(Gec,Hec*r1,Cec,Dec*r1,x0,Nm);%se mutiplica por 0 porque no hay
entrada solo con CI
V=[X0(:,4)];
U=-Ke0*X0';
r1=10;% segunda referencia
figure
stairs(k*Ts,X0)
legend('x1(k)','x2(k)','x3(k)','x4(k)')
title({'Vector de estados del sistema','CON ESTIMADOR DE ESTADOS'})
ylabel('X(kT)')
xlabel('kT(seg)')
grid
[Ye Xe]=dinitial(Gec,Hec*r1,Cec,Dec*r1,xe,Nm);
figure
stairs(k*Ts,Xe)
legend('~x1(k)','~x2(k)','~x3(k)','~x4(k)')
title({'Vector de estados estimados del sistema','CON ESTIMADOR DE ESTADOS'})
ylabel('X(kT)')
xlabel('kT(seg)')
grid
[~,E]=dinitial(Gec,Hec*r1,Cec,Dec,e0,Nm);
figure
stairs(k*Ts,E)%salida vector error de estado comparar con el integrar tiempo
establecimiento
title({'Vector de error','CON ESTIMADOR DE ESTADOS'})
grid
figure
stairs(k*Ts,V)
title({'Variable acumulacion de error','CON ESTIMADOR DE ESTADOS'})
ylabel('Variable de acumulacion del error V(kT)')
xlabel('kT(seg)')
grid
figure
stairs(k*Ts,Y0)
title({'Vector de salida','CON ESTIMADOR DE ESTADOS'})
ylabel('Y(kT)')
xlabel('kT(seg)')
grid
figure
stairs(k*Ts,U)
title({'Variable de control','CON ESTIMADOR DE ESTADOS'})
ylabel('Variable de control U(kT)')
xlabel('kT(seg)')
grid

PREGUNTA J

%%Funcion de transferencia pulso del controlador-estimador

M=ctrb(Ge0,He0);%verificación de la controlabilidad
rM=rank(M)
[Nuz, Dez]=ss2tf(Ge0,He0,Cec,Dec,1)