Control Toolbox
Control Toolbox
Control Toolbox
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)
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
roots([1 2 0 5 10])
ans =
0.8550 + 1.4809i
0.8550 - 1.4809i
-2.0000 + 0.0000i
-1.7100 + 0.0000i
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
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
damp(Gtf)
Pole Damping Frequency Time Constant
(rad/seconds) (seconds)
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)
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
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
Gret=tf([1],[4 1],'InputDelay',2)
1
exp(-2*s) * -------
4s+1
hasdelay(Gret)
1
Gret.inputdelay
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
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 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
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 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
s3
G( s)
s 2s 2
2
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
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º
e2s
Gr ( s ) 10
4s 1
G=tf([10],[4 1])
G.outputdelay=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)
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=[M.GainMargin' M.GMFrequency']
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
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á
sc
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
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%
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
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)
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
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)
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
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')
Octave
Es la versión GNU (gratuita) de Matlab y es levemente diferente, puede encontrarse en
https://www.gnu.org/software/octave/
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