TUTORIAL DE
ANÁLISIS Y CONTROL DE
SISTEMAS USANDO MATLAB
Manuel Vargas Villanueva
Este tutorial está basado en un trabajo original de:
Manuel Berenguel Soria y Teodoro Álamo Cantarero
Prepare for your exams
Study with the several resources on Docsity
Earn points to download
Earn points by helping other students or get them with a premium plan
Guidelines and tips
Prepare for your exams
Study with the several resources on Docsity
Earn points to download
Earn points by helping other students or get them with a premium plan
Tutorial de Matlab para Control
Typology: Study notes
2018/2019
1 / 44
This page cannot be seen from the preview
Don't miss anything!
Related documents
Partial preview of the text
Download Tutorial de Matlab para Control and more Study notes Rent Control in PDF only on Docsity!
TUTORIAL DE
ANÁLISIS Y CONTROL DE
SISTEMAS USANDO MATLAB
Manuel Vargas Villanueva
Este tutorial está basado en un trabajo original de:
Manuel Berenguel Soria y Teodoro Álamo Cantarero
Contenido
1 AN ´ALISIS Y CONTROL DE SISTEMAS USANDO MATLAB 1
1.1 Introducci´on................................... 1
1.2 Tratamiento mediante funciones de transferencia. Sistemas con- tinuos....................................... 1
1.2.1 Dominio Temporal............................. 2
1.2.2 Dominio Frecuencial............................ 5
1.2.3 Comandos relacionados con operaciones de bloques........... 9
1.2.4 Lugar de las ra´ıces............................. 10
1.3 Estudio temporal y frecuencial de sistemas de primer y segundo orden........................................ 12
1.3.1 Sistemas de primer orden......................... 12
1.3.2 Sistemas de segundo orden........................ 15
1.3.3 An´alisis del efecto de un cero en la respuesta temporal de un sistema de segundo orden.............................. 23
1.3.4 Influencia de polos adicionales. Polos dominantes............ 26
1.4 Tratamiento mediante funciones de transferencia. Sistemas dis- cretos....................................... 28
1.5 Tratamiento mediante descripci´on en el espacio de estados.... 30
1
Cap´ıtulo 1
AN ´ALISIS Y CONTROL DE
SISTEMAS USANDO MATLAB
1.1 Introducci´on
En lo que sigue, se va a realizar una introducci´on a los comandos de matlab relacionados con la teor´ıa de control de sistemas. Casi todas las funciones que se describen pertenecen al Con- trol System Toolbox. Las funciones principales se van a explicar sobre ejemplos demostrativos, con el fin de que su uso y comprensi´on sean lo m´as sencillos posible. Se realizar´an ejemp- los tanto en el dominio temporal, continuo y discreto, como en el frecuencial. Tambi´en se mostrar´a la forma de dar la descripci´on de un sistema lineal mediante funci´on de transfer- encia, conjunto de polos y ceros, o variables de estado. Asimismo se comentar´a la forma de manipular una funci´on de transferencia como un objeto.
1.2 Tratamiento mediante funciones de transferencia. Sis-
temas continuos
Este apartado muestra el uso de algunas de las herramientas con las que cuenta matlab para el dise˜no y an´alisis de sistemas de control. Para el ejemplo se va a partir de una descripci´on de la planta en forma de funci´on de transferencia:
H(s) =
. 2 s^2 +. 3 s + 1 (s^2 +. 4 s + 1)(s + .5)
En matlab las funciones de transferencia se introducen dando el par de polinomios numerador- denominador:
1
2 Tratamiento mediante funciones de transferencia. Sistemas continuos
num = [.2 .3 1]; den1 = [1 .4 1]; den2 = [1 .5];
El polinomio del denominador es el producto de dos t´erminos. Para obtener el polinomio resultante se usa el producto de convoluci´on (o de polinomios).
den = conv(den1,den2)
Para ver los polos (o los ceros) de la funci´on de transferencia, podemos usar: roots(den) (roots(num)). Una forma m´as completa de convertir una funci´on de transferencia dada por dos polinomios numerador y denominador, en un conjunto de factores de grado 1, correspon- dientes a los polos (z 1 , z 2 , z 3 ) y ceros (c 1 , c 2 ), de la forma:
H(s) =
K (1 − (^) cs 1 ) (1 − (^) cs 2 ) (1 − (^) zs 1 ) (1 − (^) zs 2 ) (1 − (^) zs 3 )
es mediante el comando tf2zp:
[ceros,polos,gan] = tf2zp (N,D);
que devuelve un vector conteniendo los ceros de la funci´on de transferencia, un vector con- teniendo los polos, y un escalar correspondiente a la ganancia est´atica. La funci´on comple- mentaria a ´esta tambi´en existe:
[N,D] = zp2tf (ceros,polos,gan);
Como curiosidad, cabe mencionar que el nombre de estas funciones es bastante descriptivo: ”tf-two-zp” procede de transfer-function to zero-pole form.
1.2.1 Dominio Temporal
La respuesta ante un escal´on a la entrada se puede analizar en sistemas que tengan una descripci´on en forma de funci´on de transferencia o una representaci´on en el espacio de estados, generando un vector de tiempos y usando la funci´on step:
t = [0:.3:15]’; y = step(num,den,t); plot (t,y);
4 Tratamiento mediante funciones de transferencia. Sistemas continuos
−0.2 0 5 10 15
−0.
0
0.8 Respuesta a un impulso
tiempo(seg)
Figura 1.2: Respuesta al impulso
ramp = t; y = lsim (num,den,ramp,t); plot (t,y,t,ramp); title (’Respuesta a una rampa’); xlabel (’tiempo(seg)’);
La respuesta a la rampa puede verse en la Fig. 1.3.
(^00 5 10 )
5
10
15
20
25
30 Respuesta a una rampa
tiempo(seg)
Figura 1.3: Respuesta a la rampa
Otro ejemplo caracter´ıstico es la respuesta a un ruido uniforme aleatorio. Recordamos que rand(m,n) genera una matriz m × n de n´umeros aleatorios uniformemente distribuidos entre 0 y 1. Es importante darse cuenta del filtrado que se produce en la se˜nal cuando pasa por el sistema (hace de filtro paso bajas).
noise = rand (size(t));
AN ALISIS Y CONTROL DE SISTEMAS USANDO MATLAB´ 5
y = lsim (num,den,noise,t); plot (t,y,t,noise); title (’Respuesta a un ruido aleatorio’); xlabel (’tiempo(seg)’);
La respuesta al ruido puede verse en la Fig. 1.4.
(^00 5 10 )
1
1.4 Respuesta a un ruido aleatorio
tiempo(seg)
Figura 1.4: Respuesta a ruido aleatorio
1.2.2 Dominio Frecuencial
Respuesta en frecuencia
La respuesta en frecuencia de los sistemas se puede obtener usando las funciones bode, nyquist y nichols. Si no se le ponen argumentos a la izquierda, estas funciones generan las gr´aficas por s´ı solas. En caso contrario, vuelcan los datos en los vectores de salida oportunos. A continuaci´on se presentan ejemplos de las tres posibles sintaxis de la funci´on bode:
1.- bode(num,den)
2.- [mag,phase,w] = bode (num,den)
3.- [mag,phase] = bode (num,den,w)
La primera de ellas produce un gr´afico con la magnitud en decibelios (dB) y la fase en grados. En las otras dos la magnitud se devuelve en el vector mag y est´a expresada en unidades absolutas, no en dB. Por su parte, la fase, devuelta en el vector phase, sigue siendo en grados. En esta forma, la representaci´on es m´as interactiva, en el sentido de que ”pinchando” con
AN ALISIS Y CONTROL DE SISTEMAS USANDO MATLAB´ 7
Real Axis
Imaginary Axis
Nyquist Diagrams
−2.5−1.5 −1 −0.5 0 0.5 1 1.5 2
−
−1.
−
−0.
0
1
2
2.5 From: U(1)
To: Y(1)
Figura 1.6: Diagrama de Nyquist
dB (en el eje de ordenadas) frente a fase del bucle abierto en grados (en eje de abcisas), o llamar a la funci´on sin argumento de salida. Si se quiere en forma de ´abaco, se puede usar el comando ngrid.
nichols (num,den,w), ngrid; [mag,phase] = nichols (num,den,w);
−40 −360 −270 −180 −90 0
−
−
−
0
10
20
30
40
fase del bucle abierto (grados)
Diagrama de Nichols
6 db
3 db
1 db
0.5 db
0.25 db
0 db
−1 db −3 db −6 db −12 db −20 db
−40 db
ganancia del bucle abierto (dB)
Figura 1.7: Diagrama de Nichols
M´argenes de estabilidad
Como es bien sabido en la teor´ıa cl´asica del control, los m´argenes de estabilidad son el margen de fase y el margen de ganancia. Estos m´argenes se calculan usando el comando margin.
margin (num,den);
8 Tratamiento mediante funciones de transferencia. Sistemas continuos
[mg,mf,wmg,wmf] = margin (num,den);
Como tercer par´ametro de entrada se le puede pasar el rango de frecuencias deseado. En el primer caso indicado, en el que no se le especifican par´ametros de salida, se realiza la representaci´on del diagrama de Bode, junto con una indicaci´on, mediante l´ıneas verticales de los puntos donde se mide cada uno de los m´argenes y los valores de los mismos. En la segunda variante, como salidas se obtienen los valores de los m´argenes de ganancia (no en dB), el margen de fase (en grados) y sus correspondientes frecuencias. Si existen varias frecuencias de corte marca los m´as desfavorables (ver Fig. 1.8).
−50 10 −2 10 −1 100 101
0
50
Frequency (rad/sec)
Gain dB
Gm=3.448 dB, (w= 1.287) Pm=11.64 deg. (w=1.182)
10 −2^10 −1^100
0 − − − − Frequency (rad/sec)
Phase deg
Figura 1.8: M´argenes de fase y ganancia
Efectos de los retardos
Los retardos existen en numerosas aplicaciones de control autom´atico. En sistemas lineales continuos invariantes en el tiempo, el retardo viene representado por e−sT^. La forma m´as sencilla de manipular los retardos en matlab es en el dominio de la frecuencia. N´otese que e−jwT^ = 1|−wT^. Por tanto los retardos dejan la magnitud invariable y afectan al desfase, tendiendo a inestabilizar al sistema controlado. Para prop´ositos de representaci´on mediante el comando bode, todo lo que habr´a que hacer es restar la fase del retardo a la de la funci´on de transferencia. Por ejemplo:
num = [0.2 0.3 1]; den = [1 0.9 1.2 0.5]; T = 1; % Tiempo de retardo puro w = logspace (-2,1,100)’; [mag,fase] = bode (num,den,w); faseDelay = fase - (Tw180/pi); % Sustrae la fase tras convertirla a grados subplot(211); semilogx (w, 20*log10(mag)); grid; subplot(212); semilogx (w, [fase, faseDelay]); grid;
10 Tratamiento mediante funciones de transferencia. Sistemas continuos
u
S
S
y
Figura 1.11: Conexi´on en realimentaci´on
Para operaciones de bloques m´as complejas, resulta m´as adecuado usar la herramienta simulink que tambi´en se explicar´a en estas notas.
1.2.4 Lugar de las ra´ıces
El an´alisis mediante el lugar de las ra´ıces se puede obtener definiendo un vector de ganancias deseadas (que es el par´ametro usado habitualmente para ver la evoluci´on de los polos en bucle cerrado). La sintaxis es:
r = rlocus (N,D,K); rlocus (N,D);
En la primera forma, calcula el lugar de las ra´ıces de 1 + K N D^ ((ss)) = 0, para un vector de ganan- cias especificado, K. rlocus devuelve una matriz r con length(K) filas y length(den) columnas, conteniendo la localizaci´on de las ra´ıces complejas. Cada fila de la matriz cor- responde a una ganancia del vector K. El lugar de las ra´ıces puede ser dibujado con plot(r,’x’).
En la segunda forma, que es la usada habitualmente, la funci´on directamente dibuja el lugar de las ra´ıces. Adem´as, como vemos, no es imprescindible indicar un vector de ganancias. Para la funci´on de transferencia que ven´ıamos utilizando en los ejemplos, se obtendr´ıa el lugar de las ra´ıces mostrado en la Fig. 1.
Un comando muy ´util como complemento a rlocus es rlocfind, cuya sintaxis general es:
[K,polos] = rlocfind (num,den)
Antes de introducir dicho comando es necesario haber dibujado el lugar de las ra´ıces. Al introducir rlocfind, se pide que se seleccione con el rat´on un punto determinado del lugar,
AN ALISIS Y CONTROL DE SISTEMAS USANDO MATLAB´ 11
−4 −4 −3 −2 −1 0 1 2 3 4
−
−
−
0
1
2
3
4
Eje real
Eje imaginario
Figura 1.12: Lugar de las ra´ıces
proporcionando como resultado la ganancia K en dicho punto y la localizaci´on de los polos correspondientes a esa ganancia.
Para mejorar la precisi´on, puede realizarse un zoom en torno a una zona de inter´es en la gr´afica, antes de ejecutar rlocfin. Otra utilidad para ver la sobreoscilaci´on que correspon- der´ıa a un par de polos complejos conjugados situados en el lugar, ser´ıa dibujar los lugares geom´etricos de factor de amortiguamiento (δ) y frecuencia natural (ωn) constantes, mediante el comando sgrid.
Ejemplo: Los siguientes comandos dibujan el lugar de las ra´ıces de una funci´on de trans- ferencia, realizando una conveniente ampliaci´on y resaltando las l´ıneas de factor de amor- tiguamiento 0. 5 , 0. 6 , 0 .7 y de frecuencia natural 0. 5 rad/s (Fig. 1.13).
N = 1;
D = [1 3 2 0];
rlocus(N,D); sgrid ([0.5:0.1:0.7],0.5); axis([-2.5,1,-3,3]);
Las nuevas versiones de matlab van m´as all´a y ofrecen una herramienta interactiva para ir viendo en vivo las variaciones que sufre el lugar de las ra´ıces, conforme se a˜naden, eliminan y mueven los polos y ceros de bucle abierto. A esta herramienta se accede mediante el comando rltool. Merece la pena que el alumno dedique unos minutos a esta herramienta de gran valor did´actico.
AN ALISIS Y CONTROL DE SISTEMAS USANDO MATLAB´ 13
(^00 1 2 3 4 5 6 7 8 9 )
1 Respuesta a un escalón unitario
tiempo(seg)
Figura 1.14: Respuesta a escal´on unitario del sistema de primer orden
Las dos caracter´ısticas fundamentales de un sistema de primer orden son su ganancia est´atica K y su constante de tiempo τ. La constante de tiempo es el tiempo que tarda en alcanzar el 63% de la salida. La ganancia est´atica es el cociente entre la amplitud de salida y la de entrada en el r´egimen permanente. Estos valores se pueden comprobar directamente en la gr´afica o analizando el vector de datos resultante:
yRP = ye(length(ye)); % Valor en regimen permanente n = 1; while ye(n) < 0.63yRP n=n+1; end % Constante de tiempo (0.1 es el intervalo transcurrido entre dos medidas, % se le resta 1, porque los indices empiezan en 1): tauEstim = 0.1(n-1);
La respuesta a una rampa unitaria de entrada para nuestro sistema de primer orden se puede simular mediante:
ramp = t; yr = lsim (num,den,ramp,t); plot (t,yr,t,ramp); title (’Respuesta a una rampa’); xlabel (’tiempo(seg)’); grid;
El resultado se muestra en la Fig. 1.15, en la que aparecen, tanto la rampa de entrada como la salida. El error de seguimiento en posici´on en r´egimen permanente puede obtenerse, al igual que se ha hecho anteriormente, midiendo directamente en la gr´afica o bien a partir del vector de datos resultante.
14 Estudio temporal y frecuencial de sistemas de primer y segundo orden
(^00 1 2 3 4 5 6 7 8 9 )
1
2
3
4
5
6
7
8
9
10 Respuesta a una rampa
tiempo(seg)
Figura 1.15: Respuesta a rampa unitaria del sistema de primer orden
La respuesta impulsional se puede obtener del mismo modo, pero usando en este caso la funci´on impulse (Fig. 1.16), que tiene una sintaxis similar al comando step.
yi = impulse (num,den,t); plot (t,yi); title (’Respuesta a un impulso’); xlabel (’tiempo(seg)’);
(^00 1 2 3 4 5 6 7 8 9 )
1 Respuesta a un impulso
tiempo(seg)
Figura 1.16: Respuesta a impulso del sistema de primer orden
Para finalizar con la secci´on de sistemas de primer orden, se va a comprobar que los resultados coinciden con los esperados en la teor´ıa. Como es bien sabido, los sistemas lineales de primer orden de ganancia unidad invariantes en el tiempo tienen las siguientes caracter´ısticas (nos remitimos a la bibliograf´ıa):
Respuesta a escal´on unitario: ye 2 (t) = 1 − e(−t/τ^ ), (t ≥ 0)
16 Estudio temporal y frecuencial de sistemas de primer y segundo orden
(^00 1 2 3 4 5 6 7 8 9 )
1 Respuesta teórica a escalón unitario
(^00 1 2 3 4 5 6 7 8 9 )
1
2
3
4
5
6
7
8
9
10 Respuesta te¢rica a rampa unitaria
(^00 1 2 3 4 5 6 7 8 9 )
1 Respuesta teórica a impulso
Figura 1.17: Respuesta a escal´on, rampa e impulso del sistema de primer orden
AN ALISIS Y CONTROL DE SISTEMAS USANDO MATLAB´ 17
Caso 2: Si δ = 1 → 2 ra´ıces reales iguales en SPI (l´ımite sobre-sub), sistema cr´ıticamente amortiguado.
Caso 3: Si 0 < δ < 1 → ra´ıces complejas conjugadas en SPI (subamortiguado)
Caso 4: Si δ = 0 → Respuesta oscilatoria. Sistema cr´ıticamente estable. Ra´ıces en eje imaginario.
Caso 5: Si δ < 0 → Sistema inestable, ra´ıces en SPD.
Se va a analizar el comportamiento para el conjunto de valores de δ en la respuesta a un escal´on unitario. Se supone wn = 1. Dado que el caso 3 se ver´a con m´as detalle, dado su mayor inter´es, se presentar´a en ´ultimo lugar.
Primer caso: 2 ra´ıces reales distintas (Fig. 1.18).
t = [0:0.2:20]’; wn = 1; d = 2; num = [wn^2]; den = [1,2dwn,wn^2]; ye = step (num,den,t); plot (t,ye); title (’Respuesta a un escalon unitario’); xlabel (’tiempo(seg)’); grid;
(^00 2 4 6 8 10 12 14 16 18 )
1 Respuesta a un escalón unitario
tiempo(seg)
Figura 1.18: Respuesta a escal´on del sistema de segundo orden con dos ra´ıces reales distintas
Segundo caso: 2 ra´ıces reales iguales (cr´ıticamente amortiguado) (Fig. 1.19). El sistema, en este caso el lo m´as r´apido posible, antes de hacerse subamortiguado.