República Bolivariana de Venezuela
Ministerio del Poder Popular para la Educación
Universidad Nororiental Privada Gran Mariscal de Ayacucho
Ingeniería en Sistemas
Asignatura: Sistemas de Control
Estudio de la Respuesta en Frecuencia de un
Sistema y Bucle Cerrado
Docente:
Vicenzo Mascia
Alumnos:
Gustavo Arocha C.I 28.578.750
Jesús Padrón C.I 30.223.981
Andrés Tamaronis C.I 28.579.942
Puerto Ordaz, enero 2024
Estudio de la Respuesta en Frecuencia de un Sistema y Bucle Cerrado. Hacer uso de
Matlab y Simulink.
Construir una función de transferencia sustituyendo los dígitos ABCDEFGF con los números
de su C.I., ajustados a la derecha y rellenando con cero por la izquierda si es necesario.
1AB* s+2CD
C.I.: * * G(s) = ------------------------------
AB CDE FGH s 3 +E * s2 + FG* s+9H
Por ejemplo:
109* s+234
C.I.: 09 * 345 * 678 G(s) = ------------------------------
AB CDE FGH s 3 +5 * s2 + 67 * s+98
1. Nombre y C.I. del alumno
Jesús Padrón C.I 30.223.981
2. Crear la función de transferencia G(s)
3 0 2 2 3 9 8 1 1AB* s+2CD
C.I.: G(s) = ------------------------------
A B C D E F G H s3 +E * s2 + FG* s+9H
Sustituimos:
130s + 222
G(s) = ------------------------------
s3 + 3s2 + 98s + 91
EN MATLAB
Para ello mediante el editor de texto de Matlab se ha creado la variable ‘cedula’ y entre
paréntesis y comillas simples el valor de la variable, el cual es la cedula 30.223.981. Luego
se han declarado todas las letras correspondientes a la función de transferencia de la A
hasta la H, y cada una toma un valor tipo string double correspondiente a una posición
especifica de la cadena de caracteres de ‘cedula’.
Luego de ha declarado las variables juntas como una sola, AB es igual a A multiplicada por
10 y sumándole B, esto para poder tener los valores deseados.
Finalmente se crea la función de transferencia mediante la variable Gs igualada al comando
tf(), para declarar que será una función de transferencia, y dentro de sus paréntesis los
correspondiente valores que define el enunciado, por cada s que se encuentre sola
simplemente se coloca el numero 1.
Luego mediante la ventana de comandos podemos escribir el nombre de la variable
de la función de transferencia, el cual es Gs, y obtenemos la función de trasferencia
que creamos.
3. Obtener ceros y polos del sistema G(s)
Si el sistema solo tiene polos estables continuar con los siguientes apartados.
Si el sistema resulta con polos inestables, cambiar el signo de la parte real de esos polos
inestables y obtener una nueva función de transferencia G(s) con esa modificación. Realizar
los apartados siguientes con la nueva función de transferencia.
PROCEDIMIENTO MATEMÁTICO
Obtenemos los ceros despejando s del numerador de la función transferencia
Cero: 130s + 222 = 0
130s = −222
−222
s=
130
s = - 1,7077 Único cero de la función de transferencia G(s).
Para obtener los polos es necesario factorizar:
3 2
s +3 s +98 s+91
Utilizamos el método de Ruffini y obtuvimos:
1 3 98 91 x = - 0,9474
-0,9474 -1,9446 -91
1 2,0526 96,0554 0
Con ello obtuvimos el primer polo: s = - 0,9474
Quedando ( x +0,9474 ) ( x 2 +2,0526 x +96,0554 ) de aquí podemos obtener los demás
polos restantes mediante el uso de la formula cuadrática:
−b ± √ b2−4 ac
x=
2a
Donde a = 1, b = 2,0526 y c = 96,0554
−2,0526 ± √( 2,0526 ) −4 .1(96,0554)
2
x=
2. 1
−2,0526 ± √ 4,2132−384,2216
x=
2
−2,0526 ± √−380,0084
x=
2
−2,0526 ± √ 380,0084 i
x=
2
El denominador se puede repartir para cada uno de los numeradores
−2,0526 √ 380,0084 i
x= ±
2 2
−2,0526 √ −380,0084 i
x 1= + =−1,0263+ 9,7469i
2 2
−2,0526 √−380,0084 i
x 2= − =−1,0263−9,7469 i
2 2
De aquí obtenemos los dos polos:
s=−1,0263+9,7469 i
s=−1,0263−9,7469 i
Finalmente nos queda:
𝑠 + 1,7077
𝐺(𝑠) =
(𝑠 + 0,9474)(𝑠 + 1,0263 – 9,7469𝑖)(𝑠 + 1,0263 + 9,7469𝑖)
Cero: s=−1,7077
Polos:
s=−0,9474
s=−1,0263+9,7469 i
s=−1,0263−9,7469 i
EN MATLAB
Para obtener los ceros del sistema G(s) mediante la ventana de comandos de
Matlab simplemente escribimos el comando zero(), y entre paréntesis el nombre
de nuestra variable de función de trasferencia, Gs. Para obtener los polos del
sistema G(s), escribimos el comando pole(), entre paréntesis la variable
correspondiente a nuestra función de trasferencia, Gs.
COMPARACION DE RESULTADOS
PROCEDIMIENTO
MATEMÁTICO MATLAB
CERO −1,7077 −1,7077
−1,0263+9,7469 i −1,0263+9,7469 i
POLOS −1,0263+9,7469 i −0,9474 −1,0263+9,7469 i
−0,9474
De aquí podemos observar que ambos resultados son prácticamente iguales.
Además, vemos claramente que el único cero es negativo, es decir se posiciona en
el lado izquierdo del semiplano, al igual que los polos (de los cuales 2 de ellos
poseen parte imaginaria), por ende, al encontrarse todos los polos en semiplano
izquierdo el sistema se considera estable, y tener 2 imaginarios lo hace
marginalmente estable.
4. Dibujar el diagrama de Bode y Nyquis, y obtener el valor del módulo (en
decibelios) y la fase (en grado) de la respuesta en frecuencia para w= 5 rad/seg
de G(s).
PROCEDIMIENTO MATEMÁTICO
Diagrama de Bode:
Primeramente, para graficar el diagrama de Bode debemos obtener el valor de la
ganancia estática del sistema (dominio de la frecuencia).
130 s +222
G ( s )= 3 2
s +3 s + 98 s +91
Para ello sustituimos todas las “s” por “jw” y comenzamos a resolver:
𝑠 = 𝑗𝜔
130 jω+ 222
G ( jω )= 3 2
jω +3 jω +98 jω+ 91
130 jω+222
G ( jω )= 3 2
− jω −3 jω + 98 jω+91
Luego factorizamos y agrupamos parte real y parte imaginaria quedando de la
siguiente manera:
130 jω+ 222
G ( jω )=
( 91−3 jω2 ) + jω(98− jω2)
Con ello obtuvimos la ecuación de la ganancia del sistema 𝐺(𝑗𝜔):
|G ( jω )|= √(130 ω)2 +2222
√( 91−3 ω ) +ω (98−ω )
2 2 2 2 2
Sustituimos ω por 0 para obtener el valor de la ganancia estática:
|G ( jω )|= √
2222
=2,4396
√912
Para llevarlo a dB aplicamos:
20 log 2,4396 = 7,7464 dB Aquí es donde comienza el trazo de la magnitud.
Ahora podemos comenzar a graficar los polos y ceros pasándolos al plano
imaginario.
CEROS: s1=−1,7077 → ω 1=1,7077
POLOS: s2=−0,9474 → ω 2=0,9474
s3=−1,0263 → ω 3=1,0263
s4 =−1,0263 → ω 4=1,0263
Trazamos los polos con una línea continua y los ceros con línea discontinua.
Para la Magnitud por cada polo el trazado comienza a caer 20 dB/década, primera
mente se encuentra con el polo ω 2 y luego con dos polos juntos ω 3 ω 4 , entonces
cae el doble, para medir la caída que hay de uno a otro aplicamos:
frecuencia destino
20 log
frecuencia origen
ω3 ω 4 1,0263
20 log → 20 log =0 , 69 .2=1 ,38 dB Caida de ω2 a ω3 ω 4
ω2 0,9474
Luego de ω 3 ω 4 se encuentra con un cero ω 1, se mide la caída:
ω1 1,7077
20 log → 20 log =4 , 42 dB Caida de ω3 ω 4 a ω1
ω3 ω 4 1,0263
Ahora para la Fase, tenemos que cada encuentro con un polo reduce 90º, en este
caso al encontrarse primeramente con el polo 𝒘𝟐 comienza a caer 90º.
Luego se encuentra con dos polos, por ende baja 180º.
Luego se encuentra con un cero, por lo tanto, sube 90º.
Todo esto lo podemos apreciar en el diagrama de Bode que se muestra a
continuación:
Diagrama de Nyquist:
Para poder graficar el diagrama de Nyquist se utiliza la expresión ya conocida de
Magnitud:
|G ( jω )|= √(130 ω)2 +2222
√( 91−3 ω ) +ω (98−ω )
2 2 2 2 2
Y se utiliza la forma de la fase, la cual es:
parte imaginaria de los ceros −1 parteimaginaria de los polos
∅ ( ω ) =tan−1 −tan
parte real de los ceros parte real de los polos
−1 ω ( 98−ω )
2
130 ( ω )
∅ ( ω ) =tan−1 −tan
222 91−3 ( ω2 )
Mediante esas 2 ecuaciones se sustituye cada parte correspondiente de la tabla
siguiente.
ω G(ω ) ∅ (ω )
0 2,44 0º
0.2 2,40 -5,48º
0.5 2,25 -12,11º
1 1,96 -17,43º
2 1,67 -17,69º
3 1,63 -16,17º
5 1,88 -16,34º
8 3,66 -32,42º
Ahora comenzamos a graficar comenzando con una frecuencia de 𝜔 =0, ganancia
𝐺(𝜔) = 2,44 y una fase de .
Cada una de las flechas de la gráfica representan la fase (grado de inclinación)
correspondiente a la frecuencia 𝜔 dada, y cada frecuencia 𝜔, representada con un
punto sobre cada flecha, está ubicada a una distancia = 𝐺(𝜔) correspondiente
respecto al punto 0 de los ejes (imaginario y real).
El trazo de Nyquist comienza desde el punto 𝜔 = 0 en aumento hasta conectarse e
al punto 0 de ambos ejes (imaginario y real).
Ahora para obtener el valor del módulo de la respuesta en frecuencia para
ω=5 rad / s utilizamos la ecuación de ganancia que ya obtuvimos y sustituimos ω
por 5:
|G ( jω )|= √(130 ω)2 +2222
√( 91−3 ω ) +ω (98−ω )
2 2 2 2 2
|G ( jω )|= √(130 .5)2 +2222
√( 91−3(5) ) +(5) (98−(5) )
2 2 2 2 2
|G ( jω )|=1,88001
Lo llevamos a dB:
20 log 1,88001=5,4832 dB Este es el Modulo para una frecuencia de 5 rad/seg.
Para obtener la fase de la respuesta en frecuencia para 𝝎 = 5 utilizamos la
expresión de fase ya conocida y sustituimos 𝜔 por 5 nuevamente:
−1 ω ( 98−ω )
2
−1 130 ( ω )
∅ ( ω ) =tan −tan
222 91−3 ( ω2 )
−1 5 ( 98−5 )
2
−1 130 ( 5 )
∅ ( ω ) =tan −tan
222 91−3 ( 5 2)
∅ ( ω ) =−16,3470 °
EN MATLAB
Diagrama de Bode:
Para obtener el diagrama de Bode, mediante el editor de texto de Matlab utilizamos
el comando bode(), y entre paréntesis la variable que corresponde a nuestra
función de transferencia Gs.
Aquí podemos observar como el modulo comienza en 7.71 dB.
Además dando clic sobre el trazo del Módulo y la Fase, llevándolo a una frecuencia
de 5 rad/seg tenemos un módulo de 5.52 dB y una Fase de -16.4 grados.
Diagrama de Nyquist:
Para obtener el diagrama de Nyquist, mediante el editor de texto de Matlab
utilizamos el comando nyquist(), y entre paréntesis la variable que corresponde a
nuestra función de transferencia Gs.
COMPARACION DE RESULTADOS
PROCEDIMIENTO
MATEMÁTICO MATLAB
Modulo en
respuesta a 5,4832 dB 5 , 52 dB
𝜔 = 5 𝑟𝑎𝑑/𝑠𝑒𝑔
Fase en
respuesta a −16,3470 ° −16 , 4 °
𝜔 = 5 𝑟𝑎𝑑/𝑠𝑒𝑔
Ganancia
7,7464 dB 7 , 71 dB
Estática
De aquí podemos observar que ambos resultados tienen gran similitud sin embargo
la diferencia es poca, eso demuestra que el procedimiento a mano a pesar de ser
menos preciso se acercó bastante al software Matlab. Incluso en los diagramas
Bode y de Nyquist la perspectiva general de ambos diagramas es bastante parecida
teniendo en cuenta que uno fue hecho a mano alzada y el otro mediante un software
sofisticado.
5. Obtener el valor del módulo (en decibelios) en el pico de resonancia y la
frecuencia de resonancia de G(s). (ATENCIÓN, puede que el pico de resonancia
no exista o no coincida con el valor máximo en módulo del diagrama de Bode.
Indicar si suceden estas situaciones).
PROCEDIMIENTO MATEMÁTICO
Tanto para obtener el pico de resonancia como la frecuencia de resonancia primero
debemos llevar la función de transferencia G(s) a sus formas elementales:
130 s+222
G ( s )= 2
(s+ 0,9474)(s +2,0526 s+96,0554 )
222 0,9474 96,0554 130 s +222
G ( s )=
96,0554 s+ 0,9474 s +2,0526 s+96,0554
2
222
Ahora podemos hacer uso de la ecuación general de segundo orden analizando su
representación:
2
ωn
G ( s )= 2 2
s +2 ζ ωn s +ω n
Tenemos: (s2 +2,0526 s+ 96,0554)
2
Donde ω n=96,0554
ω n=9,8008
2 ζ ω n s=2,0526 Despejamos ζ quedando :
2,0526
ζ=
2ω n
2,0526
ζ= =0,1047
2 (9,8008)
Ya que este valor se encuentra por debajo de 0,7 eso quiere decir que el
sistema al que pertenece presenta pico de resonancia
Una vez hecho esto podemos proceder a calcular lo requerido.
Frecuencia de resonancia:
La ecuación para conocer en qué frecuencia ocurre el pico de resonancia, lo
que se conoce como frecuencia de resonancia es:
ω r=ωn √ 1−2 ζ 2
ω r=9,8008 √ 1−2(0,1047)2
ω r=9,6928 rad /seg
Pico de resonancia:
Para conocer en cuál es el valor más alto alcanzado por el modulo
(magnitud), conocido como pico de resonancia, se utiliza la siguiente
expresión:
M r=20 log
| 1
2 ζ √ 1−ζ 2 |
|
M r=20 log
1
2( 0,1047) √ 1−(0,10472 ) |
M r=20 log (4,8019)
M r=13,6283 dB
EN MATLAB
Para ello simplemente obtenemos el diagrama de Bode de Gs nuevamente,
damos click derecho sobre el diagrama, presionamos la opción “Characteristics” y
luego “Peak Response”.
Eso nos devolverá esta información, donde podemos observar la frecuencia donde se da
el pico de resonancia (frecuencia de resonancia) “At Frequency (rad/s): 9,69”, y el valor
del módulo (pico de resonancia) “Peak Gain: 16,3 dB”.
COMPARACIÓN DE RESULTADOS
PROCEDIMIENTO
MATLAB
MATEMÁTICO
Frecuencia de 9,6928 rad /seg 9 , 69 rad / seg
resonancia
Pico de 13,6283 dB 16 , 3 dB
resonancia
Se puede observar que la Frecuencia de resonancia es prácticamente la misma, sin
embargo, el pico de resonancia del procedimiento matemático hay un poco de diferencia
con el calculado en Matlab.
6. Obtener la función de transferencia M(s) resultado de realimentar unitariamente y
negativamente a G(s).
M(s)
X(s) + Y(s)
K G(s)
-
Aquí tenemos que M(s) es igual a la relación entre la salida y la entrada, de esta
manera:
Y (s) K G( s)
M ( s )= =
X (s) 1+ K ( G ( s ) ) (Realimentación)
Ya que la realimentación es unitaria entonces el denominador se multiplica por 1.
K (130 s+ 222)
s +3 s2 +98 s+91
3
M ( s )=
K (130 s+ 222)
1+ 3 2
(1)
s +3 s +98 s+91
Luego resolvemos el denominador, multiplicamos la fracción por el 1 de realimentación,
quedando prácticamente igual. Después igualamos el primer 1 al denominador de la otra
fracción entre el mismo, quedando así:
K (130 s+222)
s +3 s 2+ 98 s +91
3
M ( s )=
K (130 s+ 222)
1+ 3 2
s +3 s +98 s+91
K (130 s +222)
s +3 s 2+ 98 s+ 91
3
M ( s )= 3 2
s +3 s +98 s+ 91 K (130 s+222)
3 2
+ 3 2
s +3 s +98 s+ 91 s +3 s + 98 s +91
Resolvemos esa suma de fracciones:
K (130 s+222)
s + 3 s 2+ 98 s+ 91
3
M ( s )= 3 2
s +3 s +98 s+ 91+ K (130 s +222)
3 2
s + 3 s + 98 s+ 91
Aplicamos doble C, quedando:
K (130 s+222) 3 2
s +3 s +98 s +91
M ( s )= 3 2
. 3 2
s + 3 s + 98 s+ 91 s + 3 s + 98 s+ 91+ K (130 s +222)
Cancelamos primer denominador con segundo numerador y así obtenemos la función
de transferencia del sistema M(s).
K ( 130 s +222 )
M ( s )= 3 2
s + 3 s + 98 s+ 91+ K ( 130 s+ 222 )
7. Obtener los polos del sistema G(s) en bucle cerrado para K=5
Para obtener los polos e bucle cerrado de G(s) para K=5, es necesario conocer cómo queda G(s)
en bucle cerrado:
M(s)
X(s) + Ι (s ) Y(s)
O(s)
K G(s)
-
O(s)
Tenemos que G ( s )=
Ι (s )
Pero ya que el sistema es bastante sencillo, simplemente tenemos que:
Ι ( s )=X ( s ) ,
O ( s )=Y ( s )
Por ende la función de transferencia G(s) = M(s)
K ( 130 s+222 )
G ( s )= 3 2
K=5
s +3 s + 98 s +91+ K ( 130 s +222 )
5 ( 130 s+222 )
G ( s )= 3 2
s +3 s + 98 s +91+5 ( 130 s +222 )
650 s+1110
G ( s )= 3 2
s +3 s + 98 s +91+650 s+1110
650 s+1110
G ( s )= 3 2
s +3 s + 748 s+ 1201
Ahora procedemos a factorizar mediante el Método Ruffini:
1 3 748 1201 X = -1,6105
-1,6105 -2,2377 -1201
1,3895 745,7623 0
Con ello obtuvimos el primer polo x=−1,6105
Quedando (x +1,6105)( x 2+1,3895 x +745,7623) de aquí podemos obtener los
demás polos restantes mediante el uso de la formula cuadrática:
−b ± √ b2−4 ac
x=
2a
Donde a = 1, b = 1,3895 y c = 745,7623
−1,3895 ± √( 1,3895 ) −4 . 1(745,7623)
2
x=
2.1
−1,3895 ± √ 1,9307−2983,0492
x=
2
−1,3895 ± √−2981,1185
x=
2
−1,3895 ± √2981,1185 i
x=
2
−1,3895 √ 2981,1185 i
x= ±
2 2
−1,3895 √ 2981,1185 i
x 1= + =−0,69475+27,2998 i
2 2
−1,3895 √2981,1185 i
x 2= − =−0,69475−27,2998 i
2 2
Finalmente tenemos que para G(s) en bucle cerrado con una ganancia de K=5, un total
3 polos:
s=−1,6105
s=−0,69475+27,2998 i
s=−0,69475−27,2998 i
EN MATLAB
Para realizar esto nos vamos al editor de texto de Matlab y creamos la variable K, la igualamos a 5,
luego igualamos nuestra variable de la función de transferencia Gs a
ella misma multiplicada por K, después creamos una nueva variable para la función de
transferencia realimentada, a la cual llamamos Gs2 y la igualamos al comando feedback(Gs,1),
este comando para indicar que poseer una realimentación unitaria, por eso el 1.
Una vez hecho esto nos vamos a la ventana de comandos y escribimos el comando pole(Gs2),
para obtener los polos del sistema realimentado al cual llamamos Gs2. A continuación podemos ver
el procedimiento mencionado y lo obtenido:
COMPARACIÓN DE RESULTADOS
PROCEDIMIENTO
MATEMÁTICO MATLAB
−0,69475+27,2998 i −0 , 6948+27,2998 i
−0,69475−27,2998 i −0 , 6948−27,2998i
Polos −1,6105 −1 , 6104
Mediante el procedimiento manual se han obtenido exactamente los mismos resultados que
Matlab, por ende, el procedimiento matemático es exacto.
De aquí podemos observar que solo dos polos poseen parte imaginaria, además vemos como
todos los polos son negativos, es decir, se encuentran en el semiplano izquierdo, por ende, el
sistema G(s) realimentado negativa y unitariamente y con una ganancia K=5, es marginalmente
estable.
8. Obtener los márgenes de ganancia y de fase y las frecuencias correspondientes a
esos márgenes para G(s) (En algunos casos puede ser infinito, indicarlo de ser así).
EN MATLAB
Simplemente mediante la ventana de comandos escribimos el comando margin(Gs2), ya que
nuestra función de transferencia en bucle cerrado la hemos llamado Gs2.
Esto nos arrojara la siguiente información con el diagrama de Bode,donde podemos observar en la
parte superior que el margen de ganancia es infinito (Gm=inf dB), es decir que la fase nunca corta
los -180º. Además vemos que el margen de fase es de 4,43º (Pm=4,43 deg) y se puede apreciar el
momento en el que la Ganancia corta los 0 dB el cual es a una frecuencia de 37,3 rad/seg (at 37,3
rad/s)
9. Con Simulink simular el sistema del punto 6.
MEDIANTE MATLAB CON SIMULINK
Nos vamos a la línea de comandos en matlab y creamos un array con numerador y
denominador, igualado al comando linmod(), entre paréntesis escribimos el archivo que
creamos de simulink llamado 'Modelo', damos enter, obteniendo lo siguiente:
Para observar de mejor manera la función de transferencia creamos la variable Ms, y la
igualamos al comando tf(num,den) y esto nos muestra nuestra función de
transferencia del sistema de bloques (suponiendo que la ganancia K=1).
COMPARACION DE RESULTADOS
Tanto la función de transferencia obtenida mediante el procedimiento manual como la
obtenida mediante Matlab son válidas y funcionales, a pesar de que mediante Matlab
no se aprecia la constante de la ganancia expresada como K, sino que
automáticamente supone que es 1 y arroja directamente el resultado siendo K=1.
Podemos detallar esa alteración en el denominador, donde el 228s anteriormente era
98s, y el 313 antes era 91. Esto ha sucedido debido a la realimentación negativa y
unitaria dada, además de poseer una ganancia de 1 como ya se mencionó.
CONCLUSION
Ambos procedimientos nos ha dejado en claro que tanto el sistema Gs en lazo abierto
como en lazo cerrado con K=5, son marginalmente estable al solo contar con polos
negativos y más de 1 en eje imaginario. Además, ambos procedimientos nos han dado los
valores de Fase y Ganancia de respuesta en frecuencia 5 rad/seg bastante acertados y
confiables.
Además de la función de transferencia M(s) resultado de realimentar unitaria y
negativamente a G(s), ha sido la misma obtenida en ambos procedimientos.
Tanto el procedimiento matemático como el de Matlab han arrojado resultados
considerablemente correctos, sin embargo, el procedimiento mediante Matlab al ser
realizado mediante el software sofisticado, arroja resultados mucho más precisos y de
manera más rápida.
Los resultados obtenidos mediante el procedimiento manual, al solo haberlos trabajado
considerando apenas 4 decimales, pueden ser el origen de tal diferencia respecto a
Matlab ya que este software toma en cuenta una cantidad sumamente superior a 4
decimales, por ende, es una fuente de mayor confianza.