Respuesta en el tiempo de un Sistema de Control

UNIVERSIDAD AUTÓNOMA DE NUEVO LEÓN FACULTAD DE INGENIERÍA MECANICA Y ELÉCTRICA Respuesta en el tiempo de un Sistema de Control La respuesta de un s...
1 downloads 0 Views 268KB Size
UNIVERSIDAD AUTÓNOMA DE NUEVO LEÓN

FACULTAD DE INGENIERÍA MECANICA Y ELÉCTRICA

Respuesta en el tiempo de un Sistema de Control

La respuesta de un sistema de control, o de un elemento del sistema, está formada de dos partes: la respuesta en estado estable y la respuesta transitoria. La respuesta transitoria es la parte de la respuesta de un sistema que se presenta cuando hay un cambio en la entrada y desaparece después de un breve intervalo. La respuesta en estado estable es la respuesta que permanece después de que desaparecen todos los transitorios.

Salida

t Transitorio

Estado estable

Señales de prueba típicas. Las señales de prueba que se usan regularmente son funciones escalón, rampa, parábola, impulso, senoidales, etc. Con estas señales de prueba, es posible realizar con facilidad análisis matemáticos y experimentales de sistemas de control, dado que las señales son funciones del tiempo muy simples.

INGENIERÍA DE CONTROL

1

M.C. JOSÉ MANUEL ROCHA NÚÑEZ M.C. ELIZABETH GPE. LARA HERNÁNDEZ

UNIVERSIDAD AUTÓNOMA DE NUEVO LEÓN

FACULTAD DE INGENIERÍA MECANICA Y ELÉCTRICA

Respuesta en el tiempo de un sistema de control Sea el siguiente sistema de control

La función de transferencia del sistema es C (s ) = G (s ) R (s )

C (s ) = G (s )R(s ) La respuesta en el tiempo C (t ) es obtenida tomando la transformada de Laplace inversa de C (s )

C (t ) = L

L

-1

C (s ) = L

-1

[G(s )R(s )]

-1

Respuesta en el tiempo de un sistema de primer orden

C (s ) 1 = R(s ) Ts + 1

C (s ) =

1 R (s ) Ts + 1

Respuesta al escalón unitario

La entrada escalón unitario es R (s ) =

1 s

La respuesta en el tiempo es C (t ) = L

INGENIERÍA DE CONTROL

-1

1

− t  1 1 T = 1 − e  Ts + 1 s 

2

M.C. JOSÉ MANUEL ROCHA NÚÑEZ M.C. ELIZABETH GPE. LARA HERNÁNDEZ

UNIVERSIDAD AUTÓNOMA DE NUEVO LEÓN

FACULTAD DE INGENIERÍA MECANICA Y ELÉCTRICA

Constante de tiempo, es el tiempo que tarda el sistema en alcanzar del 63.2% de su cambio total. t =T

Conforme más pequeña es la constante de tiempo la respuesta del sistema es más rápida. Tiempo de estabilización, o tiempo de respuesta es el tiempo que necesita la curva de respuesta para alcanzar la línea de 2% del valor final, o cuatro constantes de tiempo. ts = 4 T

Respuesta al impulso unitario de un sistema de primer orden

La entrada impulso unitario es

R(s ) = 1 La respuesta en el tiempo es C (t ) = L

INGENIERÍA DE CONTROL

-1

1

 1  1 −T t  Ts + 1 = T e

3

M.C. JOSÉ MANUEL ROCHA NÚÑEZ M.C. ELIZABETH GPE. LARA HERNÁNDEZ

UNIVERSIDAD AUTÓNOMA DE NUEVO LEÓN

FACULTAD DE INGENIERÍA MECANICA Y ELÉCTRICA

Forma general de la función de transferencia de primer orden C (s ) G (T1 s + 1) = R (s ) Ts + 1 Donde G es la ganancia del sistema.

Polos Son los valores de s que hacen que el polinomio del denominador sea cero. Son las raíces del polinomio del denominador. Ceros Son los valores de s que hacen que el polinomio del numerador sea cero. Son las raíces del polinomio del numerador. El Polo de la función es

s=−

1 T

s=−

1 T1

El cero de la función es

Ubicación del polo y cero del sistema en el plano s.

INGENIERÍA DE CONTROL

4

M.C. JOSÉ MANUEL ROCHA NÚÑEZ M.C. ELIZABETH GPE. LARA HERNÁNDEZ

UNIVERSIDAD AUTÓNOMA DE NUEVO LEÓN

FACULTAD DE INGENIERÍA MECANICA Y ELÉCTRICA

Ejemplo Un circuito eléctrico RC representado a continuación tiene las siguientes funciones de transferencia.

Las funciones de transferencias considerando e0 (t ), i(t ), eR (t ) como salidas y ei (t ) como entradas son: E0 (s ) 1 = Ei (s ) RCs + 1

I (s ) Cs = Ei (s ) RCs + 1

E R (s ) RCs = Ei (s ) RCs + 1

Si el valor de la resistencia es R = 100 KΩ y el valor del capacitor es C = 1 µ f y se le aplica un voltaje de entrada ei = 10V . Obtener la respuesta en el tiempo para cada salida. La constante de tiempo del circuito es: RC = 100 000 * 0.000001 = 0.1 seg

El voltaje de entrada aplicado es Ei (s ) =

10 s

El voltaje de salida E 0 (s ) sería  1  10   10  10  E0 (s ) =    =     0.1s + 1  s   s + 10  s 

Obteniendo la transformada inversa de Laplace

(

)

e0 (t ) = 10 1 − e −10 t volts

El tiempo de estabilización para una banda del 2% sería t S = 0.4 seg INGENIERÍA DE CONTROL

5

M.C. JOSÉ MANUEL ROCHA NÚÑEZ M.C. ELIZABETH GPE. LARA HERNÁNDEZ

UNIVERSIDAD AUTÓNOMA DE NUEVO LEÓN

FACULTAD DE INGENIERÍA MECANICA Y ELÉCTRICA

La corriente I (s ) sería.  1 *10 −6 s  10   10 *10 −6 s  10    =    I (s ) =      0.1s + 1  s   s + 10  s 

Obteniendo la transformada inversa de Laplace

(

)

i (t ) = 0.1 *10 −3 e −10 t ampers

El voltaje E R (s ) sería.  0.1s  10   s  10  E R (s ) =    =     0.1s + 1  s   s + 10  s 

Obteniendo la transformada inversa de Laplace

(

eR (t ) = 10 e −10 t

)

El tiempo de estabilización para una banda del 2% sería t S = 0.4 seg

INGENIERÍA DE CONTROL

6

M.C. JOSÉ MANUEL ROCHA NÚÑEZ M.C. ELIZABETH GPE. LARA HERNÁNDEZ

UNIVERSIDAD AUTÓNOMA DE NUEVO LEÓN

FACULTAD DE INGENIERÍA MECANICA Y ELÉCTRICA

Respuesta en el tiempo de un sistema de segundo orden El diagrama de bloques de un sistema de segundo orden en términos de la relación de amortiguamiento ζ y frecuencia natural no amortiguada ω n es

La función de transferencia de lazo cerrado es

ω n2 C (s ) = 2 R (s ) s + 2ζω n s + ω n2 La ecuación característica nos da información sobre el comportamiento dinámico del sistema. Las raíces de la ecuación característica s 2 + 2ζω n s + ω n2 = 0 serían

s1, 2 = −ζω n ± ω n ζ 2 − 1 Si ζ > 1 los polos de lazo cerrado son reales y diferentes, el sistema se denomina sobreamortiguado y su respuesta transitoria es exponencial s1 = −ζω n + ω n ζ 2 − 1

s 2 = −ζω n − ω n ζ 2 − 1

Si ζ = 1 los polos de lazo cerrado son reales e iguales, el sistema se denomina críticamente amortiguado y su respuesta es exponencial s1 = −ω n

s 2 = −ω n

Si 0 < ζ < 1 los polos de lazo cerrado son complejos conjugados, el sistema se denomina subamortiguado y su respuesta es oscilatoria s1 = −ζω n + jω n 1 − ζ 2

s 2 = −ζω n − jω n 1 − ζ 2

Si ζ = 0 los polos de lazo cerrado son imaginarios, el sistema se denomina sin amortiguamiento y su respuesta tiene oscilaciones mantenidas s1 = jω n INGENIERÍA DE CONTROL

s 2 = − jω n

7

M.C. JOSÉ MANUEL ROCHA NÚÑEZ M.C. ELIZABETH GPE. LARA HERNÁNDEZ

UNIVERSIDAD AUTÓNOMA DE NUEVO LEÓN

FACULTAD DE INGENIERÍA MECANICA Y ELÉCTRICA

1. Caso sobreamortiguado ζ > 1 Para una entrada escalón unitario R(s ) = 1 s

La transformada inversa de Laplace

2. Caso críticamente amortiguado ζ = 1 Para una entrada escalón unitario R(s ) = 1 s

La transformada inversa de Laplace

3. Caso subamortiguado 0 < ζ < 1 Para una entrada escalón unitario R(s ) = 1 s

ω n2

C (s ) =

(s + ζω n + jω d )(s + ζω n − jω d )s

en donde ω d = ω n 1 − ζ 2 es la frecuencia natural amortiguada

La transformada inversa de Laplace

4. Caso sin amortiguamiento ζ = 0 Para una entrada escalón unitario R(s ) = 1 s C (s ) =

ω n2

(s + jω n )(s − jω n )s

La transformada inversa de Laplace c(t ) = 1 − cos ω n t INGENIERÍA DE CONTROL

8

M.C. JOSÉ MANUEL ROCHA NÚÑEZ M.C. ELIZABETH GPE. LARA HERNÁNDEZ

UNIVERSIDAD AUTÓNOMA DE NUEVO LEÓN

FACULTAD DE INGENIERÍA MECANICA Y ELÉCTRICA

Familia de curvas c(t) con diversos valores de ζ .

Definiciones de las especificaciones de respuesta transitoria.

Con frecuencia, las características de desempeño de un sistema de control se especifican en términos de la respuesta transitoria para una entrada escalón unitario, dado que ésta es fácil de generar. (Si se conoce la respuesta a una entrada escalón, es matemáticamente posible calcular la respuesta para cualquier entrada.) Al especificar las características de la respuesta transitoria de un sistema de control para una entrada escalón unitario, es común especificar lo siguiente:

1. Tiempo de retardo, t d 2. Tiempo de levantamiento, t r 3. Tiempo pico, t p 4. Sobrepaso máximo, M p 5. Tiempo de asentamiento, t s

1. Tiempo de retardo, t d : es el tiempo requerido para que la respuesta alcance la primera vez la mitad del valor final.

INGENIERÍA DE CONTROL

9

M.C. JOSÉ MANUEL ROCHA NÚÑEZ M.C. ELIZABETH GPE. LARA HERNÁNDEZ

UNIVERSIDAD AUTÓNOMA DE NUEVO LEÓN

FACULTAD DE INGENIERÍA MECANICA Y ELÉCTRICA

2. Tiempo de levantamiento, t r : es el tiempo requerido para que la respuesta pase del 10 al 90%,del 5 al 95% o del 0 al 100% de su valor final. Para sistemas subamortiguados de segundo orden, por lo común se usa el tiempo de levantamiento de 0 a 100%. Para sistemas sobreamortiguados, suele usarse el tiempo de levantamiento de 10 a 90%. 3. Tiempo pico, t p : es el tiempo requerido para que la respuesta alcance el primer pico del sobrepaso. 4. Sobrepaso máximo (porcentaje), M p : es el valor pico máximo de la curva de respuesta, medido a partir de la unidad. Si el valor final en estado estable de la respuesta es diferente de la unidad, es común usar el porcentaje de sobrepaso máximo. La cantidad de sobrepaso máximo (en porcentaje) indica de manera directa la estabilidad relativa del sistema. M p% =

c(t p ) − c(∞ ) c(∞ )

* 100%

5. Tiempo de asentamiento, t s : es el tiempo que se requiere para que la curva de respuesta alcance un rango alrededor del valor final (por lo general, de 2%) y permanezca dentro de él.

Especificaciones de la respuesta transitoria de un sistema subamortiguado

tp =

Tiempo pico, t p Tiempo de levantamiento, t r Sobrepaso máximo porcentual, M p

tr =

π −β ωd

  −πζ

 1−ζ 2  

M p % = 100e  ts =

Tiempo de asentamiento, t s

4

σ

=

π ωd

4

ζω n

= 100e −(σ ω d )π

(banda del 2%)

La constante de tiempo del sistema es 1 ζω n

INGENIERÍA DE CONTROL

10

M.C. JOSÉ MANUEL ROCHA NÚÑEZ M.C. ELIZABETH GPE. LARA HERNÁNDEZ

UNIVERSIDAD AUTÓNOMA DE NUEVO LEÓN

FACULTAD DE INGENIERÍA MECANICA Y ELÉCTRICA

Curva de M p contra ζ .

La curva de respuesta c(t) para una entrada escalón unitario, siempre permanece dentro de un par de curvas envolventes.

La respuesta c(t ) se puede obtener por medio de las curvas envolventes El número de picos en la respuesta t sería N pi cos = s tp las magnitudes de esos picos son

(

c(t p ) = 1 + e

( )

−ζω n t p

( c(3t ) = (1 + e

c(2t p ) = 1 − e

(

) c(∞) ) ) c(∞ ) ) ) c(∞)

−ζω n 2t p

(

−ζω n 3t p

p

INGENIERÍA DE CONTROL

11

M.C. JOSÉ MANUEL ROCHA NÚÑEZ M.C. ELIZABETH GPE. LARA HERNÁNDEZ

UNIVERSIDAD AUTÓNOMA DE NUEVO LEÓN

FACULTAD DE INGENIERÍA MECANICA Y ELÉCTRICA

Ejemplo 1 Considere el sistema de la figura, en el que ζ = 0.6 y ω n = 5 rad / seg . Obtenga el tiempo de levantamiento t r , el tiempo pico t p el sobrepaso máximo M p y el tiempo de asentamiento t s cuando el sistema está sujeto a una entrada escalón unitario.

A partir de los valores dados de ζ y ω n obtenemos ω d = ω n 1 − ζ 2 = 4 , y σ = ζω n = 3 . Tiempo de levantamiento t r es

ωd  4  = tan −1   = 0.93 rad 3 σ 

β = tan −1 

tr =

π − β π − 0.93 = = 0.55 seg ωd 4

Tiempo pico t p es tp =

π π = = 0.785 seg ωd 4

Sobrepaso máximo M p es M p = e (−π σ ω d ) = e (−3π 4 ) = 0.095

M p % = 0.095 * 100 = 9.5%

Tiempo de asentamiento t s para el criterio del 2% es 4 4 t s = = = 1.33 seg σ 3

INGENIERÍA DE CONTROL

12

M.C. JOSÉ MANUEL ROCHA NÚÑEZ M.C. ELIZABETH GPE. LARA HERNÁNDEZ

UNIVERSIDAD AUTÓNOMA DE NUEVO LEÓN

FACULTAD DE INGENIERÍA MECANICA Y ELÉCTRICA

Ejemplo 2 Considere el sistema de la figura. Determine el valor de k de modo que la relación de amortiguamiento ζ sea de 0.5. Luego obtenga el tiempo de crecimiento t r , el tiempo pico t p , el sobreimpulso máximo M p y el tiempo de establecimiento t s , en la respuesta a un escalón unitario.

Solución. La función de transferencia de este sistema es: Gω n2 C (s ) = 2 R(s ) s + 2ζω n s + ω n2

C (s ) 16 = 2 R(s ) s + (0.8 + 16 k )s + 16

Comparándola con la función general de 2º orden s 2 + 2ζω n s + ω n2 = s 2 + (0.8 + 16 k )s + 16

lo que nos da que 2ζω n = (0.8 + 16 k ) y ω n2 = 16 la ganancia del sistema G = 1 La frecuencia natural no amortiguada ω n . Despejando k

k=

ω n = 4 rad / seg

2ζωn − 0.8 2(0.5)(4 ) − 0.8 = = 0.2 16 16

k = 0.2

El máximo sobrepaso M p %  −π ζ

M p % = 100 e 

1−ζ 2  

 − 0.5 π

= 100e 

1−(0.5 )2  

= 16.3%

Frecuencia natural amortiguada ω d ω d = ω n 1 − ζ 2 = 4 1 − 0.5 2 = 3.464 rad / seg

El tiempo pico t p tp =

π π = = 0.907 seg ω d 3.464

El tiempo de crecimiento (levantamiento) t r β = cos −1 (ζ ) = cos −1 (0.5) = 1.047 rad

INGENIERÍA DE CONTROL

tr =

13

π − β π − 1.047 = = 0.605 seg 3.464 ωd

M.C. JOSÉ MANUEL ROCHA NÚÑEZ M.C. ELIZABETH GPE. LARA HERNÁNDEZ

UNIVERSIDAD AUTÓNOMA DE NUEVO LEÓN

FACULTAD DE INGENIERÍA MECANICA Y ELÉCTRICA

Tiempo de estabilización (asentamiento) t s ts =

4

ζω n

=

4 = 2 seg (0.5)(4)

Podemos calcular c(∞ ) para una entrada escalón unitario utilizando el teorema del valor final   1  16   = 1 c(∞ ) = lim sC (s ) = lim s 2  s →0 s →0  s + (0.8 + 16 k )s + 16  s 

Gráfica de respuesta c(t ) para una entrada escalón unitario, utilizando la envolvente. El número de picos sería N pi cos =

ts 2 = = 2.2 t p 0.907

c(∞ ) = GR1 = 1

Los valores de la respuesta c(t ) son

(

c(t p ) = c(0.907) = 1 + e

(

= 1+ e

−1.814

) c(∞)

)(1) = 1.163

(

c( 2t p ) = c(1.814) = 1 − e

(

( )

−ζωn t p

( )

−ζωn 2t p

) c(∞ )

)

= 1 − e −3.628 (1) = 0.973

INGENIERÍA DE CONTROL

14

M.C. JOSÉ MANUEL ROCHA NÚÑEZ M.C. ELIZABETH GPE. LARA HERNÁNDEZ

UNIVERSIDAD AUTÓNOMA DE NUEVO LEÓN

FACULTAD DE INGENIERÍA MECANICA Y ELÉCTRICA

Ejemplo 3 Para el sistema de la figura, determine los valores de la ganancia K y la constante de realimentación de velocidad K h para que el sobrepaso máximo en la respuesta al escalón unitario sea 0.2 y el tiempo pico sea 1 seg. Con estos valores de K y K h , obtenga el tiempo de levantamiento y el tiempo de asentamiento. Suponga que J = 1 Kg − m y que B = 1 N − m / rad / seg . C (s ) K = 2 R(s ) Js + (B + KK h )s + K K J = K  B + KK h  2 s + s + J J  

Comparándola con la forma general del sistema de segundo orden, nos queda Gω n2 B + KK h C (s ) K = 2 2ζω n = y ω n2 = 2 J J R(s ) s + 2ζω n s + ω n la ganancia del sistema es G = 1 Como la ganancia del sistema es 1 y la entrada es un escalón de magnitud 1, la salida se estabiliza en 1, c(∞ ) = GR1 = 1 . Utilizando el teorema del valor final, se puede obtener la magnitud donde se estabiliza el sistema.   1  K   = 1 c(∞ ) = lim sC (s ) = lim s 2  s →0 s →0  Js + (B + KK h )s + K  s 

(M

Como el máximo sobrepaso es 0.2, correspondería a un 20%,  −π ζ

M p % = 100e 

1−ζ 2  

= 20%

1

ζ =

= 20%

p%

π

) 1

=

2

  M p %    ln   100  

2

+1

π2

[ln(0.2)]2

= 0.456 +1

el tiempo pico t p = 1 seg tp =

π = 1 seg ωd

ω d = π = 3.1416 rad / seg

ωd = ωn 1 − ζ 2

como ω n2 =

ωn =

K entonces J

INGENIERÍA DE CONTROL

ωd 1− ζ 2

=

π 1 − (0.456 )

2

= 3.53 rad / seg

K = Jω n2 = (1)(3.53) = 12.46 N − m 2

15

M.C. JOSÉ MANUEL ROCHA NÚÑEZ M.C. ELIZABETH GPE. LARA HERNÁNDEZ

UNIVERSIDAD AUTÓNOMA DE NUEVO LEÓN

como 2ζω n =

B + KK h J

FACULTAD DE INGENIERÍA MECANICA Y ELÉCTRICA

Kh =

2ζω n J − B 2(0.456 )(3.53)(1) − (1) = = 0.178 seg K 12.46

El tiempo de levantamiento t r

β = cos −1 (ζ ) = cos −1 (0.456 ) = 1.097 rad

tr = ts =

Tiempo de asentamiento t s

π − β π − 1.097 = = 0.65 seg ωd π

4

ζω n

=

4 = 2.485 seg (0.456)(3.53)

Gráfica de respuesta c(t ) para una entrada escalón unitario, utilizando la envolvente. t 2.485 El número de picos sería N pi cos = s = = 2.485 ≈ 2 c(∞ ) = GR1 = 1 tp 1 Los valores de la respuesta c(t ) son

(

c(t p ) = c(1) = 1 + e

(

)

( )

−ζωn t p

) c(∞)

= 1 + e −1.609 (1) = 1.2 −ζω (2 t ) c(2t p ) = c(2) = 1 − e n p c(∞ )

(

= 1− e

(

−3.219

)

)(1) = 0.96

INGENIERÍA DE CONTROL

16

M.C. JOSÉ MANUEL ROCHA NÚÑEZ M.C. ELIZABETH GPE. LARA HERNÁNDEZ

UNIVERSIDAD AUTÓNOMA DE NUEVO LEÓN

FACULTAD DE INGENIERÍA MECANICA Y ELÉCTRICA

Ejemplo 4 La figura (a) muestra un sistema vibratorio mecánico. Cuando se aplica al sistema una fuerza de 2 lb (entrada escalón), la masa oscila como se aprecia en la figura (b). Determine m, b y k del sistema a partir de esta curva de respuesta. El desplazamiento x se mide a partir de la posición de equilibrio.

Solución. La función de transferencia de este sistema es X (s ) = P(s )

X (s ) 1 = 2 P(s ) ms + bs + k

1 m b k s2 + s + m m

Comparándola con la función general de 2º orden Gω n2 b k C (s ) 1 = 2 tenemos que la ganancia G = y s 2 + 2ζω n s + ω n2 = s 2 + s + 2 R(s ) s + 2ζω n s + ω n k m m lo que nos da que 2ζω n =

b k y ω n2 = m m

Dado que la entrada es un escalón de magnitud 2, tenemos P(s ) =

2 s

sustituyendo en la función de transferencia del sistema X (s ) =

1 2   ms + bs + k  s  2

utilizando el teorema del valor final para obtener el valor en el cuál se estabiliza el sistema. de la gráfica de respuesta el valor donde se estabiliza el sistema es x(∞ ) = 0.1 x(∞ ) = lim sX (s ) = X (s ) = lim s s →0

s →0

2 2 = = 0.1 pie k s ms + bs + k

(

2

)

Por tanto k = 20 lb f / pie

Como el sistema se estabiliza en 0.1 y el sobrepaso es de 0.0095, entonces M p = 9.5% que corresponde a ζ = 0.6 .

INGENIERÍA DE CONTROL

17

M.C. JOSÉ MANUEL ROCHA NÚÑEZ M.C. ELIZABETH GPE. LARA HERNÁNDEZ

UNIVERSIDAD AUTÓNOMA DE NUEVO LEÓN

FACULTAD DE INGENIERÍA MECANICA Y ELÉCTRICA

1

ζ =

=

π2   M p %   ln    100 

2

+1

1

π2

[ln(0.95)]2

= 0.6 +1

El tiempo pico t p se obtiene mediante π π π = = 2 ωd ωn 1 − ζ 0.8ω n La respuesta experimental muestra que t p = 2 seg . Por tanto tp =

ωn =

Dado que ω n2 =

3.1416 = 1.96 rad / seg (2)(0.8)

k , obtenemos m m=

20

ω

2 n

=

20

(1.96 )2

= 5.2 slug = 166 lb

b se determina a partir de b m b = 2ζω n m = (2)(0.6 )(1.96 )(5.2) = 12.2 lb f / pie / seg 2ζω n =

INGENIERÍA DE CONTROL

18

M.C. JOSÉ MANUEL ROCHA NÚÑEZ M.C. ELIZABETH GPE. LARA HERNÁNDEZ

UNIVERSIDAD AUTÓNOMA DE NUEVO LEÓN

FACULTAD DE INGENIERÍA MECANICA Y ELÉCTRICA

Grafique la respuesta en el tiempo para el siguiente sistema de control, para una entrada escalón unitario, y determine el tiempo de levantamiento t r , el tiempo pico t p el sobrepaso máximo M p y el tiempo de asentamiento t s . R(s)

INGENIERÍA DE CONTROL

25 s (s + 2 )

+-

19

C(s)

M.C. JOSÉ MANUEL ROCHA NÚÑEZ M.C. ELIZABETH GPE. LARA HERNÁNDEZ

Suggest Documents