Matlab ® Toolbox de control

 Matlab posee una completa librería para analizar y diseñar sistemas de control lineales
invariantes en el tiempo (Control Systems Toolbox)
 Para ello posee un tipo de datos llamado sistema LTI (lineal invariante en el tiempo)
 Ese tipo de datos es un objeto y como tal admite funciones de creación, modificación,
extracción de datos, conversión de tipo de modelo, etc
 Pude usar modelos definidos por funciónes de transferencia (externo) asi como modelos
en variable de estados (internos)
 Son funciones de Matlab en el sentido que pueden opcionalmente devolver valores, es
decir usar la sintaxis
[lista de valores de salida]=función(lista de parámetros) Si devuelve escalares no son
necesarios los corchetes
 Tanto los valores de salida como los parámetros pueden ser escalares, vectores o matrices
y van separados por coma
 Puede manejar funciones de transferencia matriciales, es decir sistemas MIMO
 Se accede a la ayuda con help control o doc control (mas detallada)

Funciones de creación y conversion

Creacion de sistemas (objetos) según su representación (clase)
 Sys=tf([numerador],[denominador]) Crea sistema a través de vectores de coeficientes
numerador y denominador (forma fraccional), debe tenerse en cuenta que deben estar
completos, es decir si falta algún coeficiente debe ir un 0
 Sys=zpk([ceros],[polos],[ganancia]) Crea sistema a través de vectores de ceros, polos y
ganancia (forma factorizada), , hay que tener en cuenta el signo, para sistemas estables
los polos deben ser negativos
 Sys=ss(A,B,C,D]) Crea sistema a través de las matrices del sistema dinamico en VE (forma
matricial)
 Si se tiene la forma factorizada y se quiere usar tf() se puede usar la orden
conv([],[]) ) (producto de polinomios) la cual puede anidarse y si esta en forma
fraccional y se desea poner como polo cero se usa roots(pol) que da las raíces de
un polinomio
Ejemplo
Hallar la representación en función de transferencia polo cero y fraccional del sistema en forma
factorizada
5
G( s ) 
s( s  3)( s  2)
Gpc=zpk([],[0 -3 -2],5)
Gpc =
5
-------------
s (s+3) (s+2)
Para obtener los coeficientes del denominador se hace el producto de los factores
conv([1 0],conv([1 2],[1 3]))
1 5 6 0
Luego
Gtf=tf([5],[1 5 6 0])
Gtf =
5
-----------------
s^3 + 5 s^2 + 6 s

Ejemplo
Si se tuviera un sistema en forma racional donde falta un coeficiente debe ir 0 en su ubicacion
2s  3
G( s ) 
s  2s3  5s  10
3

Gtf=tf([2 3],[1 2 0 5 10])

Gtf =
2s+3
----------------------
s^4 + 2 s^3 + 5 s + 10
Se hallan las raices
roots([2 3])
ans =
-1.5000

roots([1 2 0 5 10])
ans =
0.8550 + 1.4809i
0.8550 - 1.4809i
-2.0000 + 0.0000i
-1.7100 + 0.0000i

Gpc=zpk([-3/2],[0.855+j*1.4809 0.855-j*1.4809 -2 -1.71],2)

Gpc =
2 (s+1.5)
------------------------------------
(s+1.71) (s+2) (s^2 - 1.71s + 2.924)

 Falta verificar que sean iguales, lo cual se hara mas adelante

Ejemplo
Hallar la representación en variables de estado del sistema
0 1 0
x(t )    x ( t )    u (t )
 1 2  5 
y(t )  1 0 x(t )  0 u (t )
A=[0 1;-1 -2], B=[0;5],C=[1 0],D=0
Gss=ss(A,B,C,D)
Gss =
a=
x1 x2
x1 0 1
x2 -1 -2
b=
u1
x1 0
x2 5
c=
x1 x2
y1 1 0
d=
u1
y1 0

Conversion de representacion
Las mismas órdenes anteriores sirven para convertir formatos de sistema poniendo como
parámetro el sistema original
 FT=tf(Sys) Convierte a función de transferencia
 ZP= zpk(Sys) Convierte a polo, cero, ganancia
 SS=ss(Sys) Convierte a Variable de estados,
 Debe tenerse en cuenta que al no haber un única representación en VE, Matlab
toma una representación llamada "companion", que en general no coincide con
las formas canónicas controlable u observable
 Smin=minreal(Sys) Realizacion minima: elimina estados no controlables ni observables en
VE o cancela polos y ceros en funciones de transferencia
Ejemplo
Verificar que sean iguales las representaciones de
2s  3
G( s ) 
s  2s 3  5s  10
4

Gpc=zpk([-3/2],[0.855+j*1.4809 0.855-j*1.4809 -2 -1.71],2)

Gtf1=tf(Gpc)
Gtf1 =
2s+3
--------------------------------------
s^4 + 2 s^3 - 1.019e-05 s^2 + 5 s + 10
2
Solo difieren en que el termino de s no es exactamente 0 pero es un valor muy pequeño 105
Ejemplo
Sea el sistema
 Y (s) 5
G( s )   2
U ( s ) s  2s  1

Las forma canonia controlable es

0 1 0
x(t )    x ( t )    u (t )
 1 2  5 
y(t )  1 0 x(t )  0 u (t )
Obtener una representacion con Matlab
Gtf=tf([5],[1 2 1])
Gtf =
5
-------------
s^2 + 2 s + 1
Gss=ss(Gtf)
Gss =
a=
x1 x2
x1 -2 -1
x2 1 0
b=
u1
x1 2
x2 0
c=
x1 x2
y1 0 2.5
d=
u1
y1 0
La matriz A esta rotada y la ganancia 5 se divide entre 2 en B y 2.5 en C, se puede verificar que son
equivalentes haciendo
Gtf1=tf(Gss)
Gtf1 =
5
-------------
s^2 + 2 s + 1

Ejemplo
Realización mínima, sea la función
s 1
G( s ) 
s  6s 2  11s  6
3

La cual en realidad tiene una cancelación polo y cero que no se ve a simple vista, y se puede
simplificar con minreal()
den=[1 6 11 6], num=[1 1]
G=tf(num,den)
zpk(G)
Gmin=minreal(G)
zpk(Gmin)

Extracción de datos
 [numerador,denominador]=tfdata(Sys) Obiene numerador y denominador de la FT de un
sistema
 [ceros,polos,ganancia]=zpkdata(Sys) Idem con polos, ceros y ganancia
 (A,B,C,D]=ssdata(Sys) Idem con las matrices en VE
 Devuelven los datos como estructura (struct) si se desea tenerlos como vectores
debe agregarse el parámetro opcional 'v'
Ejemplo
Hallar numerador y dernominador de la FT y polos, ceros y ganancia de la representación en
variables de estado del sistema
0 1 0
x(t )    x ( t )    u (t )
 1 2  5 
y(t )  1 0 x(t )  0 u (t )
A=[0 1;-1 -2], B=[0;5],C=[1 0],D=0
Gss=ss(A,B,C,D)
[n,d]=tfdata(Gss,'v') % devuelve como vectores
n=
0 0 5
d=
1 2 1
[z,p,k]=zpkdata(Gss,'v') % devuelve como vectores
z=
Empty matrix: 0-by-1
p=
-1
-1
k=
5

Funciones de análisis en el dominio temporal

Obtencion de caracteristicas
 P=pole(Sys) Devuelve los polos de un sistema
 Z=zero(Sys) Devuelve los ceros de un sistema, de no haberlos devuelve un vector vacio
  damp(Sys) Muestra los polos, el amortiguamiento (damping) y la frecuencia natural de
cada uno de ellos, además la constante de tiempo asociada
 pzmap(Sys) Dibuja el mapa de polos y ceros
 isstable(Sys) Devuelve 1 si el sistema es estable y 0 si no lo es
 Kc=dcgain(Sys) Devuelve la ganancia estatica es decir Kc  lim Sys( s )
s 0
 No confundir la ganancia estatica Kc con la del formato polo cero K , además
puede dar infinita (inf) si el sistema tiene integradores puros
Ejemplo
Obtener polos, ceros, damping, frecuencias, constantes de tiempo y ganancia estatica del sistema
con polos complejos conjugados, además dibujar el mapa de polos y ceros
5( s  1)
G( s ) 
s  8s 2  30s  36
3

Gtf=tf([5 5],[1 8 30 36])

Gtf =
5s+5
-----------------------
s^3 + 8 s^2 + 30 s + 36
pole(Gtf)
-3.0000 + 3.0000i
-3.0000 - 3.0000i
-2.0000 + 0.0000i
zero(Gtf)
-1
pzmap(Gtf)

damp(Gtf)
Pole Damping Frequency Time Constant
(rad/seconds) (seconds)

-2.00e+00 1.00e+00 2.00e+00 5.00e-01

-3.00e+00 + 3.00e+00i 7.07e-01 4.24e+00 3.33e-01
-3.00e+00 - 3.00e+00i 7.07e-01 4.24e+00 3.33e-01

dcgain(Gtf)
0.1389
 Observese que K  5 , en tanto que Kc  0.139
Ejemplo
Obtener damping, frecuencias, constantes de tiempo, estabilidad y ganancia estatica del sistema
tipo 1
2
G( s ) 
s ( s  2)
Gzp=zpk([],[0 -2],2)
Gzp =
2
-------
s (s+2)

damp(Gzp)
Pole Damping Frequency Time Constant
(rad/seconds) (seconds)

0.00e+00 -1.00e+00 0.00e+00 Inf

-2.00e+00 1.00e+00 2.00e+00 5.00e-01

isstable(Gzp)
0
dcgain(Gzp)
Inf
 Observese que la constante de tiempo asociada al polo en el origen es infinita, al
igual que la ganancia Kc y que el sistema es inestable

Graficacion dominio temporal

 impulse(Sys ) Dibuja la respuesta al impulso unitario
 step(Sys ) Dibuja la respuesta al escalón unitario
 stepinfo(Sys ) Devuelve los parámetros de respuesta al escalón unitario
 lsim=(Sys,u,t) Dibuja la respuesta a una entrada u y un tiempo t arbitrarios, dados como
vectores, ambos vectores tienen que tener el mismo tamaño
 Las ordenes de dibujo se pueden llamar con parámetros de entrada y almacenan
en vectores las salidas y y los tiempos t ,para ser dibujados con la orden plot() o
analizados mediante Matlab, la sintaxis es por ejemplo [y,t]=step(Sys), plot(t,y)
 También permiten superponer las respuestas de varios sistemas con distintos
colores por ejemplo step(Sys1,Sys2,…)
 Además es posible acotar el valor de tiempo con parámetros adicionales, como
por ejemplo agregando el tiempo final por ejemplo step(Sys1,tfinal)
 Admiten todas las opciones de formato de la orden plot()
  Haciendo click derecho -> Characteristics se pueden obtener gráficamente los
valores principales como etiquetas emergentes
Ejemplo
Obtener la respuesta al escalón en rojo del sistema ya visto y determinar sus parámetros gráfica y
analíticamente
5( s  1)
G( s ) 
s  8s 2  30s  36
3

Gtf=tf([5 5],[1 8 30 36])

step(Gtf,'r'), grid on

stepinfo(Gtf)
RiseTime: 0.2612
SettlingTime: 2.1411
SettlingMin: 0.1260
SettlingMax: 0.2022
Overshoot: 45.5997
Undershoot: 0
Peak: 0.2022
PeakTime: 0.7599
Ejemplo
Superponer las repuestas del anterior en azul (color por defecto) y el sistema con ganancia en
continua unitaria en rojo, verificar la ganancia en continua
2
G( s ) 
s  2s  2
2

Ggu=tf([2],[1 2 2])
dcgain(Ggu)
1
step(Gtf,Ggu,'r'), grid on
legend('Sistema Ganancia NO unitaria','Sistema Ganancia unitaria')

Ejemplo
Obtener y dibujar el error en la respuesta a la rampa del sistema con ganancia en continua unitaria
2
G( s ) 
s  2s  2
2

t1=0:0.01:10; % genera tiempo

u=t1; % entrada 0 tiempo (rampa)
[y,t]=lsim(Ggu,u,t1); % simula y no dibuja
e=u-y' % error
ess=e(end);
1.0000
plot(t,e) , grid on
title('Error a la rampa')
Retardos de tiempo puros y aproximaciones racionales
 Sys = tf(num,den,'InputDelay',retardo) Define una función de transferencia con retardo a
la entrada
 Sys = tf(num,den,'OutputDelay',retardo) Idem a la salida
 [Link]=retardo Agrega o modifica el retardo a la salida
 [Link]=retardo Idem a la entrada
 [Link] Devuelve el retardo a la salida
 [Link] Idem a la entrada
 hasdelay(Sys) Devuelve 1 si tiene retardo y 0 si no lo tiene
 [num,den] = pade(T,N) Devuelve el numerador y denominador de la aproximación de
racional de Pade de orden N del retardo T

 En sistemas SISO es indistinto el retardo a la entrada o a la salida

Ejemplo
Crear un sistema con retardo de la forma
e2s
Gr ( s ) 
4s  1
Dibujar su respuesta y comparala con la aproximación de pade de orden 1

Gret=tf([1],[4 1],'InputDelay',2)
1
exp(-2*s) * -------
4s+1
hasdelay(Gret)
1
[Link]
2
[n,d]=pade(2,1)
Gpade=tf(n,d)
-s + 1
------
s+1
step(Gret)
hold on
G=tf([1],[4 1])
step(G*Gpade,'r')
step(Gpade,'g')
legend('FT con retardo puro','FT con aproximacion de Pade','Aproximacion de Pade')
Conexión de sistemas
 Sys = series(Sys1,Sys2,….SysN) o Sys = Sys1*Sys2*…. *SysN) Conexión en serie
 Sys = parallel(Sys1,Sys2,… SysN) o Sys = Sys1+Sys2+… +SysN) Conexión en paralelo

 Ambas pueden incluir constantes

 F = feedback(G,H) Conexión de lazo cerrado por realimentacion negativa de la rama

G( s )
directa G con la rama realimentación H, equivale a F ( s ) 
1  G( s ) H ( s )
 Para realimentacion positiva se hace F = feedback(G,-H) equivale a
G( s )
F ( s) 
1  G( s ) H ( s )
 Sysc = connect(sys1,...,sysN,operación,entradas,salidas) Conexión general (puede
buscarse en la documentación), generalmente con las anteriores basta y no será tratada
aquí, operación puede definirse con la orden operación=sumblk(formula)
Ejermplo
Hallar Y (s) / R(s) del sistema de lazos anidados

Se definen las FT
G1=tf([1],[1 1])
1
-----
s+1
G2=tf([1],[1 0])
1
--
s
Se resuelve el lazo interno

G3=feedback(G1,1)
1
-----
s+2

El producto es
G4=10*G3*G2
10
---------
s^2 + 2 s
F=feedback(G4,1)
F=
10
--------------
s^2 + 2 s + 10

Ejemplo
Hallar Y (s) / R(s) del sistema con bloques en paralelo, verificar estabilidad y hallar la ganancia
estatica

Se definen las FT y constantes

Kp=5
Cd=zpk([0],[],2)
G=zpk([],[0 -10],1)
H=3

Se resuelve el paralelo
C=Kp+Cd
Se resuelve la serie y la realimentación
F=feedback(C*G,H)
2 (s+2.5)
-------------------
(s+1) (s+15)
isstable(F)
1
dcgain(F)
0.3333

Lugar de raices
 rlocus(Sys) Dibuja el Lugar de Raíces de un sistema K*Sys para K positiva en una ventana
interactiva que permite desplazarse sobre el mismo, hacer zoom, cambiar el trazo, etc.
 rlocus(Sys,k) Idem con un vector de ganancias k definido por el usuario
 rlocus( Sys1,Sys2….,SysN) Dibuja el Lugar de Raíces de varios sistemas (pueden usarse
colores y tipos de línea definidos para cada uno como en plot())
 K = rlocfind(Sys,P) Obtiene el valor de ganancia K para la ubicación de un polo en P
 [K,Polos] = rlocfind(Sys) Superpone un cursor sobre el lugar de raíces y devuelve las
ganancias y los polos correspondientes
 Aquí la sentencia grid o sgrid produce una grilla con lugares de amortiguamiento y
frecuencia natural constantes
 Si se desea LR con ganancia negativa se hace rlocus(-Sys)
Ejemplo
Dibujar el LR de la función sin y con grilla y hallar la ganancia y ubicación del polo doble
2 s 2  5s  1
G( s)  K 2
s  2s  3
G = tf([2 5 1],[1 2 3])
rlocus(G)
sgrid % dibuja grilla

[K,polos]=rlocfind(G) % seleccionar el punto de separacion

K=
0.9604
Polos=
-1.1644 + 0.0091i
-1.1644 - 0.0091i
Verificar los polos de lazo cerrado
F=feedback(K*G,1), pole(F)
-1.1644 + 0.0091i
-1.1644 - 0.0091i
Ejemplo
Dibujar el LR de la función y hallar ganancia y puntos de separación y Ganancia Critica de
estabilidad
1
G( s)  K
s( s  2)( s  3)
G=zpk([],[0 -2 -3],1), rlocus(G)

Aquí la ganancia critica de estabilidad Kcr=29 y la separación se produce con Ks=2.11 con un polo
doble en p=-0.785 y otro polo real menor a -3
Para hallar la ubicación del polo restante se hace
[K,P]=rlocfind(G,-0.785)
K=
2.1126
P=
-3.4305
-0.7850
-0.7845

Es decir el tercer polo para ganancia 2.11 se encuentra en -3.43 y su relación es

P(1)/P(2)
4.3701
Por lo que habría una leve dominancia
Ejemplo
Mediante el LR se puede diseñar un controldor proprcional, supóngase que el el ejemplo anterior
se desea ubicar los polos para que el Sobrepaso sea de 5%
G=zpk([],[0 -2 -3],1), rlocus(G), sgrid % moverse hasta encontrar ubicacion

Es necesaria una ganancia de 3.64 y los polos dominantes están ubicados en -0.69+j 0.721
Luego el tiempo de establecimiento estará dado por
K=3.64, P1=-0.69+j*0.721
tss=4/(-real(P1))
5.7971
El tiempo del primer pico se puede calcular como
tp=pi/imag(P1)
4.3573
El polo restante esta en
[k,p]=rlocfind(G,P1)
k=
3.6111
p=
-3.6173 + 0.0000i
-0.6914 + 0.7213i
-0.6914 - 0.7213i
Es decir será real y estará aproximadamente en -3.6173, se vera la relación de partes reales
real(p(1))/real(p(2))
5.2321
Por lo que puede ser considerado dominante, se verificara ahoa la función de lazo cerrado
F=feedback(K*G,1), step(F), grid on

Puede verse que los valores no difieren mucho de los esperados pero la diferencia se debe al
tercer polo real
Ejemplo
Contorno de raíces, cuando en el sistema varia un parámetro aparte de la Ganancia se pueden
superponer lugares de raíces para valores discretos de ese parámetro mediante un lazo for, sea la
función con un polo variable en -a
1
G(s)  K , a  1, 2,3
( s  5)( s  a)
n=1
for a=[ 1 2 3 ]
switch n
case 1
spec='b' % a=1 azul
case 2
spec='r' % a=2 rojo
case 3
spec='k' % a=3 negro

end
G=zpk([],[0 -5 -a],1)
rlocus(G,spec), hold on
n=n+1
end

Puede verse que la Kcr varia, para a=1 es 34.6, para a=2 es 70.2 y para a=3 es 119, es decir cuanto
mas alejado del origen esta el polo, mayor es la ganancia critica

Ejemplo
Dibujar el LR para ganancia negativa (en rojo) y positiva (en negro línea discontinua) en un solo
grafico de la función y hallar la ganancia Ks y ubicación del polo doble y la ganancia critica Kcr de
existir
s3
G( s) 
s  2s  2
2

G = tf([1 3],[1 2 2])

rlocus(-G,'r',G,'--k')
Se puede ver que los LR son complementarios para ganancia negativa la Ks=0.472 con polo doble
en P=-0.774 y la Kcr=0.665, para ganancia positiva la Ks=8.47 con polo doble en P=-5.24 siendo
incondicionalmente estable
Identificacion
Supongase que se tiene tabulados los valores de respuesta de un sistema a un escalon unitario, en
este caso en un archivo . mat que contiene las variables t e y
Se desea conocer el orden y los parámetros del sistema, esto se llama identificación

Sistema de segundo orden monótono con polo dominante

Es de la forma
KC
G( s) 
1s  1 2 s  1
Ejemplo
load ident_1 t y
Aparecen las variables en el workspace
plot(t,y)
Evidentemente es de 2do orden pues es en forma de "S" y KC  0.9 entonces el valor de tiempo
de establecimiento es y(tss )  0.9  0.98  0.88 lo cual se puede buscar con el cursor como se
muestra dando tss  2.1
Tambien se puede buscar por código usando indices haciendo
Kc=y(end) % Ganancia continua
erel=(Kc-y); % error relativo
its=find(abs(erel)<0.02,1,'first') % busca el índice donde el error relativo se hace <2%
84
tss=t(its) % saca el valor de tiempo de establecimiento
 1.9111
erel(its) % verifica el valor
0.0199

taud=tss/4 % valor de tau dominante

0.4778
El sistema aproximado se puede poner entonces como
Ga=tf([Kc],[taud 1])
0.898
----------------
0.4778 s + 1
Y se comparan ambas respuestas no encontrando diferencias significativas, solo al inicio
[ya,ta]=step(Ga); plot(t,y,ta,ya,'r'), legend('Resp medida', 'Resp aproximada'),grid

Sistema de segundo orden Subamortiguado

Ejemplo
load ident_2 t y
Aparecen las variables en el workspace
plot(t,y)

Es un sistema subamortiguado en el cual es difícil determinar gráficamente el tss que ronda los 3
seg, el tiempo de pico cerca de 2.2 seg y el máximo 1.255, se hara por codigo
yfin=y(end) % salida final
1.2001
ymax=max(y) %valor maximo
1.2551
S=(ymax-yfin)/ymax % sobrepaso
0.0439
ip=find(y==ymax) % indice del pico
34
tp=t(ip) % tiempo de pico
2.1792
wa=pi/tp % frecuencia amortiguqada
1.4416
zita=-log(S)/sqrt(pi^2+log(S)^2) % amortiguamiento
0.7054
Kc=yfin
erel=flipud((Kc-y)); % error relativo invertido en el tiempo
its=find(abs(erel)>0.02,1,'first') % busca el índice donde el error relativo se hace > 2%
61
tf=flipud(t) % tiempo invertido
tss=tf(its) % saca el valor de tiempo de establecimiento
3.03
w0=wa/sqrt(1-zita^2) % frecuencia natural
2.0339
4/(zita*w0) % debe ser parecido al tiempo de establecimiento obtenido antes
2.7878
Ga=tf([Kc*w0^2],[1 2*zita*w0 w0^2])
4.964
--------------------
s^2 + 2.87 s + 4.137
[ya,ta]=step(Ga); plot(t,y,ta,ya,'r'), legend('Resp medida', 'Resp aproximada'),grid

Son prácticamente iguales

Funciones de análisis en el dominio frecuencial

Diagramas de Bode (cartesiano)
 bode(Sys) Dibuja el diagrama de Bode de un sistema
 bode(Sys1,...SysN) Idem de varios sistemas
 bode(...,w) Idem con un vectorde frecuencias creado por el usuario
 [magnitud,fase,wout] = bode(sys) Devuelve vectores de magnitud, fase y frecuencia

Diagramas de Nyquist (polar)

 nyquist(Sys) Dibuja el diagrama de Nyquist
 nyquist(Sys1,...SysN) Idem de varios sistemas
 nyquist(….,w) Idem con un vector de frecuencias creado por el usuario
 [re,im,w] = nyquist(sys) Devuelve vectores de parte real, imaginaria y frecuencia
 Valores de la respuesta y margenes estabilidad
 BW = bandwidth(Sys) Ancho de banda -3dB del sistema
 BW = bandwidth(Sys,dbdrop) Idem con frecuencia de corte definida por el usuario
 [Gm,Pm,Wgm,Wpm] = margin(sys) Margenes de ganancia (en veces, no en dB), de fase y
frecuencias correspondientes
 margin(sys) Margenes de ganancia y de fase de forma grafica en el Bode
 S = allmargin(sys) Devuelve un objeto struct con margenes de ganancia, de fase, retardo
y frecuencias correspondientes, muy útil para sistemas con retardo
Ejemplo
Hallar el Bode, el Nyquist , obtener los márgenes de fase y ganancia del sistema en forma grafica y
en forma de datos y además la ganancia estatica
5
G( s ) 
s( s  3)( s  2)
G=zpk([],[0 -3 -2],5)
figure(1)
bode(G,'r'), grid
figure(2)
nyquist(G,'r')
[MG,MF,Wmg,Wmf] = margin(G)
MG =
6
MF =
55.1540
Wmg =
2.4495
Wmf =
0.7559
Da un MG= 6 el cual se puede convertir a dB con Matlab
MGdB=20*log10(MG)
MGdB =
15.5630
Kc=dcgain(G)
Inf
Ejemplo
Hallar el diagrama de polos y ceros, el Bode, los márgenes de ganancia y fase, la ganancia estatica
en dB, el pico de resonancia y el ancho de banda -3dB y -6 dB del sistema
s4
G( s )  1000
s  14s  145s 2  1252s  2020
4 3

G=tf(1000*[1 4],[1 14 145 1252 2020])

figure(1)
pzmap(G)
figure(2)
bode(G), grid
KcdB=20*log10(dcgain(G))
5.9342
magnitud=bode(G);
PicodB=20*log10(max(magnitud))
11.4695
BW3 = bandwidth(G)
3.4405
BW6 = bandwidth(G,-6)
 12.6687

Se observa que se cumplen los valores calculados, el ancho de banda se toma después del pico en
este caso
Ejemplo
Hallar el diagrama de Nyquist de la FT del ejemplo anterior, hallar la frecuencia y la amplitud de
los cruces por -180º y los márgenes de ganancia y estabilidad

Ejemplo
Para un sistema con retardo de la forma que se muestra, graficar el Nyquist y hallar el primer
cruce por 180º
e2s
Gr ( s )  10
4s  1
G=tf([10],[4 1])
[Link]=2
nyquist(G)
Se observa que a escala real la parte de alta frecuencia se achica
Se puede aumentar la resolucion creando un vector de ganancias personalizado desde la
frecuencias 0.12 hasta 12 con 1000 valores logaritmicamente espaciados
w=1.2*logspace(-1,1,1000);
nyquist(G,w)

Se pueden calcular todos los márgenes mediante

M=allmargin(G)

M=
GainMargin: [1x9 double]
GMFrequency: [1x9 double]
 PhaseMargin: 170.7712
PMFrequency: 2.4868
DelayMargin: 1.1985
DMFrequency: 2.4868
Stable: 0
Para crear los pares de vectores márgenes de ganancia y frecuencia se apilan en una matriz losl
campos GMFrequency y GainMargin transpuestos

MGF=[[Link]' [Link]']
0.3807 0.9183
1.5866 3.9585
2.8362 7.0862
4.0900 10.2219
5.3454 13.3611
6.6011 16.5009
7.8572 19.6413
9.1134 22.7820
10.3677 25.9181

Funciones de diseño
Sisotool
Es una herramienta interactiva de diseño que permite ajustar los polos, ceros y ganancias de un
compensador y moverlos con el mouse y/o ajustarlos escribiendo sus valores, la sintxis es
 sisotool() Abre un diseño en blanco
 sisotool(Sys) Abre un diseño tomado como proceso a Sys
Se ejecuta y abre 2 ventanas, la de trabajo que contiene el LR y el Bode, y la principal que tiene las
opciones, esta nunca se debe cerrar ya que intenta cerrar la herramienta, al hacerlo pide si se
desea guardar la tarea para después continuar, ello se guarda en Design History
sisotool()
Si se llamó sin argumentos se deben importar la Planta G, el sensor H y el prefiltro F, los cuales
deben estar definidos previamente, por defectos toma a todos como 1
En la ventana System Data se selecciona el elemento a importar, hacer click en browse y abre la
ventana Model Import que permite elegir un modelo
Diseño en el dominio temporal
Se ilustrara mediante un ejemplo, sea el proceso de tercer orden y tipo 1 (con integrador puro)
280
G( s) 
s  23,3s 2  92, 6s
3

Las especificaciones son

 Sobrepasamiento:
S  5%,
 Tiempo de asentamiento al 2%:
te2%  1seg
 Error de régimen Nulo ante entrada escalon
Se comienza definiendo la Planta en Matlab y llamando la librería con ella como argumento
G=tf([280],[1 23.3 92.6 0])
G=
280
-----------------------
s^3 + 23.3 s^2 + 92.6 s
sisotool(G)

Por comodidad se quitara el Bode haciendo View -> Design configuration y poniendo None en el
grafico del Bode
Ahora se agregan las especificaciones (Design Constraints) verificando que el tiempo de
establecimiento que se desea sea para el 2 % y que muestre el controlador en forma de polo cero
y ganancia, Matlab permite elegir estos parámetros en Edit-> Toolbox preferences

Las especificaciones se agregan haciendo Click derecho -> Design requirements ->new
Ahora la ventana muestra en color la zona que cumple los requerimientos

Los polos deseados deben estar en la intersección de las líneas, como el lugar de raíces no pasa
por ese punto un controlador P no bastara
A continuación se elegirá y diseñara el controlador que será un control de adelanto ya que no hace
falta acción integral al ser el sistema tipo 1
El controlador será
sc
C (s)  K
s p
Y de modo que el cero cancele el polo dominante de la planta, para ello se hace
pole(G)
0
-18.2168
-5.0832
Luego se deberá ubicar el cero fijo en -5.0832, ello se hace con Click derecho -> Add Pole/Zero ->
Real Zero o con el icono del cero O , se toma el cero y se lo coloca cerca del polo, luego se ajusta
con Click derecho -> Edit Compensator y en la otra ventana se selecciona y escribe el valor dando
Enter

A continuación resta agregar el polo móvil con el icono X y moverlo para que el LR pase por el
lugar deseado

Por ultimo se ajusta la ganancia (cuadrado magenta) moviéndola hasta que este sobre el lugar
deseado como se muestra
Se puede previsualizar la respuesta y la acción de control en la ventana de control haciendo en la
pestaña Analisis Plot y chequear lo siguiente, r es la entrada, y es la salida y u la señal de control
Con Show Analysis Plot se abre el LTI viewer que muestra en azul la salida y en verde la acción de
control

Por ultimo se exporta el controlador a Matlab con File -> Export

Se selecciona Compensator C y con Export to Workspace se pasa a Matlab, si existiera la variable

pide sobreescibir o cambiando Export As se le pone otro nombre
C
2.3962 (s+5.083)
------------------------
(s+10.9)
Por ultimo se realiza en Matlab la verificación de la función de lazo cerrado haciendo
F=feedback(C*G,1)
step(F)

Se aprecia que el sobrepaso es levemente menor (4.42) y el tiempo de establecimiento levemente

mayor (1.11), la razón se encuentra en el tercer polo el cual se muestra con
pzmap(F)

Diseño en el dominio frecuencial

Cuando hay sistemas con retardo el LR no funciona bien ya que aproxima por Pade el retardo, aquí
es mas conveniente diseñar en el dominio de la frecuencia sea el sistema con retardo de la forma
e2s
Gr ( s )  10
4s  1
Haciendo
G=tf([10],[4 1])
[Link]=2
sisotool(G)

Las especificaciones son

 Sobrepasamiento:
S  5%,
 La respuesta más rápida que se pueda encontrar con el error acotado a 5%
Luego el sobrepaso corresponde a  0.707 y conlleva a un margen de fase de
2
MF  t g 1
1   4  2 2
psita=0.707
MF=atand(2*psita/(sqrt(sqrt(1+4*psita^4)-2*psita^2)))
65.5246
Entonces como antes se quita esta vez el LR y se deja el Bode
Se observa que los márgenes son negativos por ende el sistema es inestable a lazo cerrado
Se planteara un controlador de atraso
sc
C (s)  K
s p
Como antes se puede poner un cero que cancele el polo de la planta
pole(G)
-0.2500

También se agrega un polo real negativo 10 veces menor, p=-0.025

Se tiene la siguiente situación
Se ajusta la ganancia hasta llegar a un MF=65º, ello se hace haciendo drag sobe el bode y
levantando o bajando

Con Menú-> Analysis -> Response to Step Command se puede visualizar la salida, la cual es

Tiene un error de 10%, subiendo la ganancia y alejando hacia la izquierda el polo se puede reducir
el error manteniendo el MF y viendo la respuesta al escalon, quedando al final
Y la respuesta al escalon tiene un error de 4%

Luego se exporta el controlador a Matlab según se vio resultando

C=
0.096853 (s+0.25)
--------------------------
(s+0.009556)
El cual resulta un compensador de atraso ya que el polo esta a la izquierda del cero y se muestran
los Bodes de ambos en azul y rojo
bode(G,C,'r')

Se observa que aporta -67 º a 0.04 rad/seg, se vera el sistema de lazo cerrado
F=feedback(C*G,1); % lo muestra en VE
step(F)

bode(F)
Sistemas Digitales
Muchas de las funciones sirven para sistemas discretos, solo hay que especificar en las ordenes de
creación el periodo de muestreo Ts como parámetro
 Sys=tf([numerador],[denominador],Ts) Crea sistema a través de vectores de coeficientes
numerador y denominador (forma fraccional), debe tenerse en cuenta que deben estar
completos, es decir si falta algún coeficiente debe ir un 0
 Sys=zpk([ceros],[polos],[ganancia], Ts) Crea sistema a través de vectores de ceros, polos
y ganancia (forma factorizada), , hay que tener en cuenta el signo, para sistemas estables
los polos deben ser negativos
 Sys=ss(A,B,C,D],Ts) Crea sistema a través de las matrices del sistema dinamico en VE
(forma matricial)
 A Partir de aquí Matlab ya lo trata como un objeto sistema LTI de tempo discreto

Diferencias
En el lugar de raíces aparece la grilla con los valores de amortiguamiento y frecuencia natural en el
plano z con la orden grid o zgrid
2
En el Bode se muestra la respuesta hasta la mitad de la frecuencia de muestreo s 
Ts
En las ordenes de dibujo de respuesta temporal se muestra la respuesta en formato retenida de
orden cero, si se desea verla en forma de barras se puede usar la orden stem
Si el sistema tiene retardo puro, es conveniente que el tiempo de muestreo sea un múltiplo entero
del mismo y no mas de 2 veces el retardo
Funciones de conversión
 Sysd=c2d(Sysc,Ts, 'metodo') Convierte de continuo a discreto
 Sysc=c2d(Sysd,Ts, 'metodo') Convierte de discreto a continuo
 Sysd1=d2d(Sysd,Ts, 'metodo') Convierte de discreto a discreto con otro periodo de
muestreo
El parámetro opcional método indica de que forma se convierte, los mas usados son
 'zoh' Con retentor de orden 0 a la entrada, es el método por default y el mas usado
 'foh' Idem de orden 1
1  sTs / 2
 'tustin' Aproximacion Bilineal z 
1  sTs / 2
 pT cT
 'matching' Mapea los polos y ceros en z e i
,z e i
pi ci
 Los 2 primeros sirven para discretizar procesos continuos y los 2 ultimos para
controladores diseñados en dominio analogico

Digitalización de Procesos continuos

En general la entrada al proceso debe ser analógica y por ello se utilizan los retentores de orden 0
y algunas veces el de orden 1

Ejemplo
Hallar el sistema discretizado y dibujar los polos y ceros con Ts  0.1 del sistema
5
G( s ) 
s( s  3)( s  2)
Con los métodos zoh y foh

G=zpk([],[0 -3 -2],5)
5
-------------
s (s+3) (s+2)
Gzoh=c2d(G,0.1)
0.00073665 (z+3.3) (z+0.236)
-------------------------------------------
(z-1) (z-0.8187) (z-0.7408)

Sample time: 0.1 seconds

Puede verse que los polos se mapean en z  e sTs y aparecen ceros que no estaban
Gfoh =
0.00018876 (z+8.979) (z+0.9049) (z+0.09118)
------------------------------------------------------------------
(z-1) (z-0.8187) (z-0.7408)
Sample time: 0.1 seconds
Aquí aparece un cero mas y los polos se mantienen
Ejemplo
Hallar el sistema discretizado, la gancia estatica y dibujar las respuestas al escalon con Ts  0.1
del sistema
5
G( s ) 
( s  3)( s  2)
Con zoh, foh
G=zpk([],[-3 -2],5), dcgain(G)
5
-----------------
(s+3) (s+2)
0.8333
Gzoh=c2d(G,0.1), dcgain(Gzoh)
0.021203 (z+0.8465)
------------------------------
(z-0.8187) (z-0.7408)
0.8333
Gfoh=c2d(G,0.1,'foh'), dcgain(Gfoh)
0.0073665 (z+3.3) (z+0.236)
----------------------------------------
(z-0.8187) (z-0.7408)
0.8333
Puede verse que la ganancia estatica y los polos no varian
Ejemplo
Para el sistema anterior dibujar las respuestas al escalon, los diagramas de Bode (en Hz) y los
lugares de raices del sistema continuo y del discretizado con Retentor de orden 0 y Ts  0.1
step(G,Gzoh,'r'), grid on, legend('Sistema Continuo', 'Sistema discreto')

Puede verse que coinciden en los instantes de muestreo

bode(G,Gzoh,'r'), grid on, legend('Sistema Continu o', 'Sistema discreto')

Puede observarse que a baja frecuencia coinciden y se alejan conforme se aproxima a la mitad de
la frecuencia de muestreo, en este caso 5 Hz
Ademas la fase del sistema discreto cruza por -180º a una frecuencia de1.6 Hz dando un margen
de ganancia de MG=26.8 dB, indicando estabilidad condicional a lazo cerrado relativa, el continuo
es incondicionalmente estable

Para poner la frecuencia en Hz click derecho -> properties

rlocus(G,Gzoh,'r'), grid on, legend('Sistema Continuo', 'Sistema discreto')

Mientras el sistema continuo es incondicionalmente estable, el discreto tiene una Ganancia critica
(cruce con el circulo unitario) de Kcr=22.9
Ejemplo
Comparar los lugares de raíces discretos del ejemplo anterior con zoh y foh
rlocus(Gzoh,Gfoh)

El foh tiene una ganancia critica Kcr=245 y el zoh es el menos estable con Kcr=21.7
Ello es debido a los ceros extra que introducen el segundo
Ejemplo
Comparar los Bodes del ejemplo anterior del sistema continuo, y con foh
bode(Gfoh,G)

Se observa que el foh aproxima mejor al sistema continuo en frecuencia que el zoh

Sisotool discreta
Se usa de la misma manera que en tiempo continuo, solo cambia en el LR que trabaja en el plano z
y el Bode es hasta la mitad de la frecuencia de muestreo
Ejemplo
Para el sistema
5
G( s ) 
( s  3)( s  2)

Digitalizar con Retentor de orden cero y Ts  0.01

Diseñar un controlador con las especificaciones

 Sobrepasamiento:
S  5%,
 Tiempo de asentamiento al 2%:
te2%  2seg
 Error de régimen Nulo ante entrada escalon
G=zpk([],[-3 -2],5)
T=0.01
Gd=c2d(G,T)
Gd =
0.00024587 (z+0.9835)
---------------------------------
(z-0.9802) (z-0.9704)
sisotool(Gd)
Lanza la herramienta y como antes se quita el bode quedando

La cual esta fuera de escala, hay que hacer zoom en la zona del 1 y entrando las especificaciones
como antes se tiene

A continuación se elegirá y diseñara el controlador que será un control PID ya que hace falta
acción integral al ser el sistema tipo 0
El controlador será de forma que un cero cancele el polo no dominante de la planta
( z  c)( z  0.9704)
CD ( z )  K
z ( z  1)
Para ello agrega ese cero y los 2 polos fijos en z=0 y z=1

A continuación se mueve el cero restante y se ajusta la ganancia para ubicar las raíces

Se exporta el controlador como Cd resultando

Cd
42.808 (z-0.9704) (z-0.9594)
----------------------------------------
z (z-1)
Esta en forma polo cero, se puede llevar a Ganancias con
Cpid=pid(Cd)
Ts 1
Kp + Ki * ------ + Kd * --------------------
z-1 Tf+Ts/(z-1)

with Kp = 2.95, Ki = 5.15, Kd = 0.399, Tf = 0.01, Ts = 0.01

Se hace el sistema de lazo cerrado y se comprueba la respuesta

F=feedback(Cpid*Gd,1); step(F)

El sobrepaso resulta levemente mayor de lo esperado y el tiempo de establecimiento un poco

menor, ello es debido a la presencia del cero del controlador circulado
Digitalizacion controladores
Si se tiene un controlador diseñado en tiempo continuo y se desea pasarlo a discreto se hace con
c2d() y se usan las opciones
 pT cT
 'matching' Mapea los polos y ceros en z e i
,z e i
es ideal para
pi ci
controladores en el dominio del tiempo
1  sTs / 2
 'tustin' Hace la transformación bilineal z  y es conveniente para el diseño en
1  sTs / 2
el dominio frecuencial

Ejemplo
Se diseño mediante Sisitool para el proceso y especificaciones
280
G( s) 
s  23,3s 2  92, 6s
3

 Sobrepasamiento:
S  5%,
 Tiempo de asentamiento al 2%:
te2%  1seg
 Error de régimen Nulo ante entrada escalon
Se obtuvo el controlador
2.4( s  5)
C ( s) 
( s  11)
Se define el proceso
G=tf([280],[1 23.3 92.6 0])
G=
280
-----------------------
s^3 + 23.3 s^2 + 92.6 s
Luego el controlador
C=zpk([-5],[-11],2.4)
2.4 (s+5)
--------------
(s+11)
Se discretizara el proceso con zoh y el controlador con matching para T=0.01
T=0.01
0.0100
Gd=c2d(G,T)
 4.405e-05 z^2 + 0.0001663 z + 3.921e-05
----------------------------------------------------------
z^3 - 2.784 z^2 + 2.576 z - 0.7922
Ahora el controlador
Cd=c2d(C,T,'matching')
2.33 (z-0.9512)
----------------------
(z-0.8958)
Los sistemas de lazo cerrado y sus respuestas son
F=feedback(C*G,1);
Fd=feedback(Cd*Gd,1);
step(F,Fd,'r')

Las respuestas son prácticamente indistinguibles

Octave
Es la versión GNU (gratuita) de Matlab y es levemente diferente, puede encontrarse en
[Link]
La mayoría de los comandos funcionan , a veces con modificaciones y no tiene (al momento) la
herramienta sisotool y el comando rlocfind entre otros
Para usar las herramientas de control hay que instalar la librería (Pakage) control
pkg Install –forge control
Y cargarla en cada sesión con pkg load control