Ejercicio Resuelto.
Sistemas de Control Automático
Ejercicio Resuelto
Jorge Pomares Baeza
Grupo de Innovación Educativa en Automática
© 2011 GITE – IEA
-1-
Índice 1. 2. 3. 4. 5. 6. 6.1 6.2 6.3 7. 8. 9. 9.1 9.2 9.1 10. 11.
4
Descripción del sistema y ecuaciones del mismo. ............................................. ............................................. 6 Diagrama de bloques y función de transferencia. ............................................ ............................................ 10 Errores y Controlador P. ................................................................ .................................................................................. .................................................. 15 Diseño de sistemas PD ..................................................................................... ..................................................... 26 Diseño de sistemas PI ................................ Diseño ................................................................ ................................ ...................................................... ...................... 34 Diseño de un controlador PID empleando la respuesta en frecuencia. ............ 38 Calcular el error de velocidad (determinar K a ) ................................ .............................................. .............. 39 Diseño de la red de avance de fase ................................................................ .................................................................. ................................ 44 Diseño de la red de retraso de fase................................................................ ................................................................... ................................ 49 Diseño de un PID por las reglas de Zigher Zigher-Nicols Nicols .......................................... .......................................... 56 Construcción de un PID discreto. ................................ ................................................................ ................................ .................................... .... 58 Análisis de los sistemas obtenidos. ................................................................ .................................................................. ................................ 77 PD ................................................................ ................................................................................................ .................................................................................... .................................................... 78 PID (Diseñado por respuesta en frecuencia).................................................... .................................................... 79 PID Discreto (Diseñado por lugar de las raíces)................................ .............................................. .............. 84 ANEXO A Valores comerciales de resistencias ................................ .............................................. .............. 87 ANEXO B Valores comerciales de condensadores ......................................... ......................................... 88
5
1. Descripción del sistema y ecuaciones del mismo.
Las partes principales de que consta todo sistema de control en bucle cerrado son la entrada (señal de referencia que índica de algún modo lo que se quiere obtener), salida (resultado o actuación en el medio que se obtiene de ejecutar el bucle de control) y ppor or último la señal de error error.. Esta señal se obtiene restando a la señal de referencia, la señal de salida de la planta, o medición del cambio que nuestro sistema ha inducido en el ella. Esta también puede estar multiplicada por un factor de corrección corrección.. Y de eesta sta manera la señal de error represent representa, a, cuanto le queda por corregir al bucle de control para obtener el resultado deseado en la salida. En el caso de nuestro sistema de control para el servomecanismo, servomecanismo, tenemos los siguientes valores de tensiones tensiones, que re representarán presentarán cada una de estas partes VR (t ) , VS (t ) y
Va (t ) respectivamente. Los dos primeros valores de tensión representan el ángulo de giro de potenciómetros potenciómetros,, mientras que el segundo es la diferencia entre los dos primeros. Enn el caso de nuestro sistema la salida ha de ser la misma que la señal de referencia, por lo que no se aplica ningún factor de corrección a la lectura de la señal de salida salida.. A esto se denomina, sistema con realimentación unitaria. Lo que si haremos es amplificar amplificar la señal de error en un factor K a . E Esto sto amplificará el valor de error aumentando así la velocidad del sistema sistema.. De esta manera, tendremos que ajustar este valor para obtener el resultado que deseamos en nuestro controlador controlador. Las ecuaciones que representan las señales comentadas serían los siguientes: VR (t ) = KPθR (t ) → Potencial de referencia ((orientación orientación deseada de la antena). VS (t ) = K Pθ S (t ) → Potencial de salida (orientación actual de la antena). Va (t ) = Ka (VS (t ) − VR (t )) = Ka K P (θS (t ) − θ R (t )) → Potencial de error. Del mismo modo, modo, la potencia de error, podemos hallarla del otro lado de la malla como la suma de la caída de tensión en la resistencia, la caída en la bobina y la corriente contra electromotriz que se produce en la bobina del motor motor.. Esta corriente es la que se genera en la bobina de un motor como re reacción acción al giro de un bobinado en el seno del campo magnético del mismo.
δθ m (t ) → Corriente contra electromotriz. δt δ i (t ) δ i (t ) δθ (t ) Va (t ) = Ra ia (t ) + La a + Vb (t ) = Ra ia (t ) + La a + K b m δt δt δt Vb (t ) = K b
6
El sumatorio de pares de fuerzas es igual al momento de inercia por las aceleraciones angulares. De este modo la ecuación que describe los pares y fuerzas que intervienen en el primer eje se puede describir como sigue a continuación:
Tm (t ) = Kmia (t ) → Par motor. TN1 (t ) y TN2 (t ) → Par generado por los engranajes que se opone al movimiento. Ti ∑ Sumatorio de los pares que intervienen
δ 2θ m (t ) = τm t 2 δ momento de inercia aceleración del eje angular
δθ (t ) δ 2θ m (t ) Tm (t ) − bm m − TN1 = τ m δ t δt2 par motor par generado par provocado por el rozamiento
TN 2 (t ) =
por la oposición al movimiento de los engranajes
− N2 TN (t ) → Relación entre los pares debido a los engranajes. N1 1
Tm (t ) = τ m
δ 2θ m (t ) δθ (t ) N + bm m + − 1 TN (t ) 2 δt δt N2 2
De igual modo en el segundo eje tenemos la siguiente ecuación:
δθ (t ) δ 2θ S (t ) TN2 (t ) − bc S = τc δ t δ t 2 Par transmitido al eje
Par que se provoca debido al rozamiento y que se opone al movimiento
Momento de inercia del segundo eje por la aceleración angular del mismo
Despejamos en la ecuación el par de entrada al eje:
δ 2θ S (t ) δθ (t ) TN (t ) = τ c + bc S 2 δt δt 2
De aquí podemos obtener la ecuación que relaciona el par del motor inicial con todo el sistema de engranajes engranajes,, y combinando ambas ecuaciones tenemos tenemos: Tm (t ) = τ m
δ 2θ m (t ) δθ m (t ) N1 δ 2θ S (t ) δθ (t ) + b + − + bc S τ c m 2 2 δt δt δt δ t N 2
A continuación vamos a simplificar esta ecuación para dejara en función función de los valores del eje motor.
7
De la relación de transmisión de los engranajes podemos obtener la equivalencia de giro del motor de salida en función del motor interno como expresamos a continuación:
θ S (t ) = −
N1 θ m (t ) N2
Dado que la derivada de una función por una constante es igual a la constante por la derivada de la función tenemos: Tm (t ) = τ m
N1 δθ m (t ) δ 2θ m (t ) δθ m (t ) N1 N1 δ 2θ m (t ) + b + − τ − + b c m c − 2 δ t2 δt N 2 δ t N 2 N 2 δ t
Tm (t ) = τ m
δ 2θ m (t ) δθ m (t ) N12 δ 2θ m (t ) N12 δθ m (t ) + b + τ + b m c c δt2 δt N22 δ t 2 N22 δ t
N 2 δ 2θ m (t ) N12 δθ m (t ) Tm (t ) = τ m + τ c 1 2 + b + b m c 2 N2 δ t 2 N2 δ t τeq ⇒ Momento de inercia equivalente
Tm (t ) = τ eq
Beq ⇒ Rozamiento equivalente (fricción viscosa)
δ 2θ m (t ) δθ (t ) + Beq m 2 δt δt
Con lo que en resumen tenemos las siguientes ecuaciones que describen el sistema: 1) V (t ) = K K [θ (t ) − θ (t ) ] a a P S R δ ia (t ) δ i (t ) δθ (t ) + Vb (t ) = Ra ia (t ) + La a + K b m 2) Va (t ) = Ra ia (t ) + La δt δt δt 3) Tm (t ) = K mia (t ) 2 4) Tm (t ) = τ eq δ θ m (t ) + Beq δθ m (t ) δ t2 δt N 5) θ S (t ) = − 1 θ m (t ) N2
8
Ahora hemos de pasar las ecuaciones al dominio de L Laplace aplace, aplace, donde una integración se traduce en una división por ‘s’ y una derivada en una multiplicación por ‘s’: ‘s’
1) V ( s) = K K [θ ( s) − θ ( s)] a a P S R 2) Va ( s) = Ra I a ( s) + La sI a ( s) + Kb sθ m ( s) = I a ( s) [ Ra + sLa ] + Kb sθ m ( s) 3) Tm ( s) = K m I a ( s) 2 4) Tm ( s) = τ eq s θ m ( s) + Beq sθ m ( s) N 5) θ S ( s) = − 1 θ m ( s) N2
9
2. Diagrama de bloques y función de transferencia. Este es el diagrama de bloques que corresponde con el sistema que estamos estudiando. θ R (s)
+
−θ
−Ka K p S (s)
Va (s) +
−
τ m ( s) Ia (s ) 1 θ m ( s) 1 Km 2 Ra + La s τ eq s + beq s
−
N1 N2
−
N1 N2
θ S (s)
θ m ( s)
Kb s
Ahora hemos de reducir el sistema sistema:: θ R (s)
+
−θ
−Ka K p S (s)
Va (s) +
Km ( Ra + La s )(τ eq s 2 + beq s )
−
θ m ( s)
θ m ( s)
Kb s
Dado un sistema en bucle cerrado como el siguiente: siguiente:
R( s)
+
−
G ( s)
Y (s) Y ( s ) = G ( S )[ R ( S ) − H ( s )Y ( S )]
H (s)
Es fácil realizar la siguiente demostración:
Y ( s ) = G ( S )[ R ( S ) − H ( s )Y ( S )] Y ( S ) = G ( s ) R ( s ) − H ( s )Y ( s )G ( s ) Y ( s ) + H ( s )Y ( s )G ( s ) = G ( s ) R ( s ) Y ( s ) [1 + H ( s )G ( s ) ] = G ( s ) R ( s ) Y (s) G (s) = R( s ) 1 + H ( s )G ( s )
10
θ S (s)
Por lo que Y ( s ) =
G( s) R( s) y con esto se obtiene la relación: 1 + H ( s)G ( s )
R( s)
+
G ( s)
−
Y (s)
R( s )
G (s) Y (s) 1 + G (s) H (s)
H (s)
De esta manera podemos podemos reducir los bucles del sistema siguiendo esta misma regla:
θ R (s)
+
−θ
−Ka K p S
−θ
−Ka K p S
−θ
Va (s)
θ m ( s)
Km ( K s) 1+ ( Ra + La s ) (τ eq s 2 + beq s ) b
K m ( Ra + La s ) (τ eq s 2 + beq s )
( Ra + La s ) (τ eq s 2 + beq s ) + K m K b s ( Ra + La s )(τ eq s 2 + beq s )
−
θ m ( s)
θ S (s)
N1 N2
−
N1 N2
θ S ( s)
( s)
θ R (s)
+
Va (s)
(s)
θ R (s)
+
Km ( Ra + La s )(τ eq s 2 + beq s )
−Ka K p S
Va (s)
Km ( Ra + La s ) (τ eq s 2 + beq s ) + K m Kb s
θ m ( s)
−
N1 N2
θ S ( s)
( s)
Así que el diagrama de bloques que representa al sistema que estamos analizando es el siguiente:
θ R (s)
+
−θ
K a K p K m N1
( Ra + La s ) (τ eq s 2 + beq s ) + K m K b s N 2 S
θ S (s)
(s)
Si aplicamos los valores de los parámetros parámetros, obtenem obtenemos os la función de transferencia del sistema sistema. 11
Los os parámetros que hemos de utilizar son los siguientes: · Ra = 2 K Ω · La = 5 H
bm = 0 N ·m / (rad / s )
·τ m = 0.05 Kg·m 2
N1 = 1
·τ c = 10 Kg·m 2
N 2 = 10
·bc = 5 N ·m / (rad / s )
K p = 24 / 2π V / rad
· K m = 0.2863 N ·m / A · K b = 0.4 V / (rad / s )
Y el resultado es el sistema: 24 ·0.2863 2π 10·( 2000 + 5s ) (τ eq s 2 + beq s ) + 0.2863·0.4s
θ R (s)
+
Ka·
−θ
S
θ S (s)
(s)
Para terminar de calcular su función de transferencia hemos de calcular el valor τ eq de y de Beq :
N12 1 τ eq = τ m + 2 τ c = 0.05 + 2 10 = 0.05 + 0.1 = 0.15 Kg m 2 N2 10 2 N 1 5 Beq = bm + bc 1 2 = 0 + 5 2 = = 0.05 N m/(rad/s) N2 10 100
θ R (s)
+
24 ·0.2863 2π 10·( 2000 + 5s ) ( 0.15s 2 + 0.05s ) + 0.2863·0.4 s Ka·
−θ
S
θ S (s)
(s)
Por lo que tenemos que la función de transferencia de nuestro sistema es:
24 ·0.2863 2π Gb ( s ) = 10·( 2000 + 5s ) ( 0.15s 2 + 0.05s ) + 0.2863·0.4 s Ka·
12
Ahora simplificaremos la ecuación para obtener una representación más sencilla del sistema, tomaremos por convenio redondear a 6 decimales:
Gb ( s) =
1.0935854 K a 10·300s 2 + 100s + 0.75s 3 + 0.25s 2 + 0.11452s
Gb ( s ) =
1.0935854 K a 3000 s 2 + 1000 s + 7.5s 3 + 2.5s 2 + 1.1452 s
Gb ( s ) =
1.0935854 K a 1001.1452 s + 3002.5s 2 + 7.5s 3
Gb ( s) =
1.0935854 K a s (1001.1452 + 3002.5s + 7.5s 2 )
Para la mayoría de los cálculos que hemos de realizar en el análisis del sistema, es aconsejable tener la función de transferencia en forma de producto de factores. En este caso particular, el numerador solo tiene término independiente por lo que hemos de descomponer en raíces el denominador. Por lo tanto lo primero que tenemos que hacer es sacar factor común en el denominador, y luego descomponer la ecuación de segundo grado resultante en un producto de raíces raíces.. Así que comenzamos normalizando el denominador: denominador
Gb ( s) = Gb ( s) =
1.0935854 K a 7.5s (133.486027 + 400.333333s + s 2 )
Gb ( s) =
13
1.0935854 K a s (1001.1452 + 3002.5s + 7.5s 2 )
0.145811K a s (133.486027 + 400.333333s + s 2 )
Ahora que ya tenemos el denominador normalizado, podemos hallar sus raíces, resolviendo la ecuación de segundo grado obtenida. obtenida
s=
−400.333333 ± (400.333333) 2 − 4 ⋅133.486027 2
s=
−400.333333 ± 160266.777511 − 533.944108 2 s=
s=
−400.333333 ± 159732.833403 2
−400.333333 ± 399.665902 −0,3337155 = 2 −399.999618
Por lo que nuestro sistema queda factorizado de esta forma:
Gb ( s ) =
14
0.145811K a s ( s + 0.337154 )( s + 399.999618 )
3. Errores y Controlador P P. Ahora vamos a analizar el sistema que hemos obtenido, desde el punto de vista de los errores y la estabilidad, para ver hasta qué punto podemos modificar el sistema para obtener los resultados que dese deseamos. amos.
G0 ( s) =
1.0935854 K a s (1001.1452 + 3002.5s + 7.5s 2 )
Los tipos de errores en un sistema de control se dividen entre los errores debidos debidos a las perturbaciones y los errores de respuesta del sistema. Entre los errores debidos a las perturbaciones tenemos; las perturbaciones de las entradas por interferencias interferencias,, lo que nos lleva a obtener salidas no deseadas, y las perturbaciones en realimentación dadas en los sensores que se emplean para realimentar realimentar sistemas en bucle cerrado, cerrado, estas últimas suelen ser debidas sobre todo al desgaste de los sistemas sistemas. Entre los errores intrínsecos de los sistemas tenemos los errores en régimen permanente como el error en velocidad y el error en aceleración, y lo loss errores en régimen estacionario como el error de posición. El sistema que estamos estudiando es un sistema en bucle cerrado, se elige este tipo de sistemas ya que son sistemas que permiten el ajuste para corregir este tipo de errores. Generalmente esta capacidad de los sistemas en bucle cerrado depende directamente del tipo de sistema del que estemos hablando. El tipo de un sistema se corresponde con el 1 número de integraciones en la función de transferencia Gb ( s ) . s En la tabla siguiente tenemos un resumen de la capacidad frente a la mitigación de errores de los sistemas en bucle cerrado, dependiendo del tipo del sistema que estemos tratando.
Tipo \
Error
Posición Velocidad
Aceleración
0
1 1+ K p
∞
∞
1
0
1 Kv
∞
2
0
0
1 Ka
El tipo de un sistema corresponde con la cantidad de raíces en cero de su denominador en la función Gb ( s ) . Por tanto een este caso tenemos un sistema de tipo 1, enn esta clase de sistemas encontramos por sus características que carecen de error de posición, tienen un error en velocidad constante y un error de aceleración infinito. Esto quiere decir que mientras que el error en posición es cero, el error en velocidad 15
podemos eliminarlo o mitigarlo mediante las acciones de control pertinen pertinentes tes,, mientras que el error de aceleración no se puede eliminar del todo en este tipo de sistemas. Los errores de un sistema se interpretan como la diferencia de la señal de salida o lectura de la planta con respecto a la señal de referencia introducida. Los distintos errores comentados anteriormente corresponden a esta diferencia en instantes muy, muy lejanos en el tiempo, pero habiendo introducido en cada uno de ellos una señal de entrada diferente. En el primer caso, el error de posición, se corresponde con eell error obtenido cuando el tiempo tiende tiende a infinito habiendo introducido una señal escalón unitario (en en el tiempo equivale a un valor de 1 constante constante) y que en el dominio de Laplace se corresponde con 1/s. En el caso del error en velocidad, se emplea una entrada en rampa para comprobar la velocidad con que el sistema responde y de esa manera comprobar el valor al que converge el error de este tipo. Por último en el caso del error en aceleración, se evalúa mediante la entrada de una parábola. Por lo tanto tanto para obtener el valor de los errores comentados hay que resolver el siguiente límite e = lim e(t ) sin embargo al intentar resolver este límite aparece una t →∞
convolución como se puede ver en la siguiente demostración:
e(t ) = r (t ) − y(t ) = r (t ) − e(t )* g (t ) Al no poder resolver el valor de e(t) de una forma sencilla por depender de una convolución, lo mejor es recurrir al teorema del valor final según el cual: lim f (t ) = lim s f ( s ) t →∞
s→0
Así pasamos a res resolver olver el problema en el dominio de Laplace donde donde la solución es más sencilla que en el dominio del tiempo. Podemos definir el error de un sistema de control con realimentación unitaria como la diferencia entre la salida y el valor de referencia E ( s ) = R ( s ) − Y ( s ) y de aquí podemos realizar la ssiguiente iguiente deducción:
E(s) = R(S ) − Y (s)
R(s) E (s)
+
G(s)
−
E(s) = R(s) − E(s)G(s)
Y ( s)
E(s) − E(s)G(S ) = R(s) E(s) (1 + G(s) ) = R(s) E ( s) =
1 R(s) 1 + G(s)
Por lo tanto según el teorema del valor final el límite que hemos de resolver para hallar el valor del error es el siguiente:
ess = lim e(t ) = lim s E(s) = lim s t →∞
16
s →0
s →0
1 R(s) 1 + G(s)
Comenzaremos analizado el error en posición del sistema, este error se define como el error en la salida respecto a la señal de referencia en un instante muy, muy lejano lejano,, siendo la entrada de referencia un escalón unitario unitario.. Sigui Siguiendo endo con esto, en el dominio de L Laplace aplace la entrada en escalón unitario se corresponde con
1 . Por tanto, el límite que s
hemos de resolver es el siguiente:
ess = lim s
1 1 1.0935854Ka s 1+ 2 s (1001.1452 + 3002.5s + 7.5s )
ess = lim s
1 1 1.0935854Ka s 11+ + 2 s (1001.1452 + 3002.5s + 7.5s )
s →0
s →0
ess = lim s →0
1 s (1001.1452 + 3002.5s + 7.5s 2 ) + 1.0935854Ka s (1001.1452 + 3002.5s + 7.5s 2 )
ess = lim s →0
s (1001.1452 + 3002.5s + 7.5s 2 )
s (1001.1452 + 3002.5s + 7.5s ) + 1.0935854Ka 2
=
0 1.0935854Ka 1.0935854K
ess = 0 Tenemos así, un error de posición igual a cero en nuestro sistema cosa que será común en todo sistema de tipo 1.
17
Respecto al error en velocidad se ha comentado que corresponderá con una constante, para llegar al valor de este error, hemos de analizar el error del si sistema stema en el tiempo al introducir una señal de referencia de tipo rampa, lo que en el dominio de 1 Laplace se corresponde con 2 . Por este motivo, el límite que hemos de resolver es el s siguiente:
ev = lim s
1 1 1.0935854Ka s2 1+ s (1001.1452 + 3002.5s + 7.5s 2 )
ev = lim s
1 1 1.0935854 K a s 21 11+ + s (1001.1452 + 3002.5s + 7.5s 2 )
s →0
s →0
1 1 s →0 1.0935854 K a s 1+ 2 s (1001.1452 + 3002.5s + 7.5s )
ev = lim
ev = lim ss→ →0
1 1 2 s (1001.1452 + 3002.5s + 7.5s ) + 1.0935854 K a s s (1001.1452 + 3002.5s + 7.5s 2 )
s (1001.1452 + 3002.5s + 7.5s 2 )
1 1001.1452 = s→0 s 1001.1452 + 3002.5s + 7.5s 2 + 1.0935854 K 1.0935854K Ka s 1.0935854 ( ) a
ev = lim
ev =
915,470525 Ka
El resultado final que tenemos es que el error en velocidad depende de K a . C Cómo ómo ya veremos más adelante tendremos que ajustar este valor para conseguir que el error en velocidad cumpla los requisitos exigidos a nuestro sistema de control control.
18
Por último analizaremos el error de aceleración para comprobar que no podemos eliminar este error como se ha comentado anteriormente. Para comprobar el valor de este error, hemos de someter al sistema a la entrada de una parábola parábola,, lo que en el tiempo 1 equivale a t 2 y en el dominio de Laplace aplace a 3 . Por lo tanto analizar esto en el infinito s equivale por el teorema del valor final a resolver el siguiente límite límite:
ea = lim s s →0
ea = lim s s →0
1 1 1.0935854Ka s3 1+ s (1001.1452 + 3002.5s + 7.5s 2 )
1 1 1.0935854Ka s32 1+ s (1001.1452 + 3002.5s + 7.5s 2 )
1 1 s →0 1.0935854Ka s2 1+ s (1001.1452 + 3002.5s + 7.5s 2 )
ea = lim
ea = lim s →0
1 1 2 s (1001.1452 + 3002.5s + 7.5s 2 ) + 1.0935854Ka s s (1001.1452 + 3002.5s + 7.5s 2 )
s (1001.1452 + 3002.5s + 7.5s 2 )
1 s →0 s 1001.1452 + 3002.5s + 7.5s 2 + 1.0935854 K s 2 1 ( ) a
ea = lim ea = lim s 0 s→
1001.1452 + 3002.5s + 7.5s 2
s (1001.1452 + 3002.5s + 7.5s 2 ) + 1.0935854Ka s
=
1001.1452 0
ea = ∞ Como hemos comprobado se cumple lo correspondiente a los sistemas de tipo 1 para cada una de las clases de error que hemos analizado analizado. alizado. Ahora hemos de estudiar el sistema desde el punto de vista de la estabilidad. Pese a que los sistemas en bucle cerrado nos permiten corregir errores, existe la posibilidad de que el sistema se convierta en inestable a causa de la realimentación. Unn sistema dinámico se considera estable, si para cualquier entrada acotada, genera como respuesta una salida también acotada. Para conocer si un sistema es estable tenemos la condición de estabilidad que nos dice que si todos los polos de la función de tr transferencia ansferencia G( s) tienen parte real negativa, el sistema es estable. Para comprobarlo tendríamos tendríamos que encontrar los polos del sistema en bucle cerrado, cerrado, que se corresponden con la lass raíces de la ecuación del polinomio característico de los sistemas en bucle cerrado 1 + G ( s) H ( s ) . Esta ecuación 19
es equivalente a 1 + G( s) en un sistema con realimentación unitaria como el nuestro nuestro.. Para evitar esto podemos usar el criterio criterio de Routh Routh-Hurwirtz Hurwirtz que nos permite determinar el número de raíces con parte real positiva de un polinomio sin la necesidad de resolverlo, por lo que podemos aplicarlo en nuestro caso caso.
20
Para aplicar el criterio de Routh Routh--Hurwirtz Hurwirtz tenemos que crear la tab tabla la de Routh RouthHurwirtz, y de ella sabremos que la cantidad de polos con parte real positiva es igual a la cantidad de cambios de signo en la primera columna de esta tabla. Para crear la tabla de RouthRouth-Hurwirtz Hurwirtz tenemos que encontrar la ecuación característica de nuestro sistema sistema,, sin embargo es más sencillo obtener la ecuación del sistema completo en bucle cerrado y comprobar los polos de esta ecuación ecuación: Gbc ( s ) =
1 1 + Gb ( s )
1.0935854 K a s (1001.1452 + 3002.5s + 7.5s 2 ) Gbc ( s ) = 1.0935854 K a 1+ s (1001.1452 + 3002.5s + 7.5s 2 ) 1.0935854 K a
Gbc ( s )
s (1001.1452 + 3002.5s + 7.5s 2 )
s (1001.1452 + 3002.5s + 7.5s 2 ) + 1.0935854 K a s (1001.1452 + 3002.5s + 7.5s 2 )
Gbc ( s )
1.093585 K a s (1001.1452 + 3002.5s + 7.5s 2 ) + 1.0935854 K a
Gbc ( s )
21
1.093585 K a 1.0935854 K a + 1001.1452 s + 3002.5s 2 + 7.5s 3
Por lo que el polinomio característico de nuestro sistema es equivalente a D( x) = 1.0935854 K a + 1001.1452s + 3002.5s 2 + 7.5s3 . Un Unaa vez que tenemos el polinomio característico podemos comenzar construyendo la tabla de Routh-Hurwirtz Routh Hurwirtz teniendo en cuenta lo siguiente: Suponemos que el polinomio característico del sistema es:
an s n + an−1s n−1 + an−2 s n−2 + ... + +a2 s 2 + a1s1 + a0 = 0 Tenemos que construir la tabla de la siguiente forma:
sn
an
an − 2
an − 4 …
an −1
an − 3
an − 5 …
s n − 2 bn −1 s n −3 cn −1
bn −3
bn −5 … siendo cn −5 …
s
n −1
⋮
cn −3
⋮
s1 g n −1 s 0 hn −1
⋮
⋮
g n −3
bn−i = −
1 an an−i −1 an−1 an−1 an−i −2
cn−i = −
1 an−1 an−i −1 bn−1 bn−1 bn−i −2
y así sucesivamente...
Además el método de Routh Routh--Hurwirtz Hurwirtz tiene la propiedad de que los coeficientes de cualquier fila se pueden multiplicar o dividir por un número positivo sin alterar los cambios en signo de la primera columna, por lo que podemos hallar el rango de valores de K a para los que nuestro sistema es estable como veremos a continuación. Primero construimos la tabla de Routh-Hurwirtz Hurwirtz a partir de nuestra función característica D( x) = 1.0935854 K a + 1001.1452s + 3002.5s 2 + 7.5s3
1001.1452 s 3 7.5 s 2 3002.5 1.0935854 K a 0 s1 b1 0 c1 0 s b1 = −
7.5 1001.1452 1 3002.5 3002.5 1.0935854 K a
b1 = −
1 ( 7.5 ⋅1.0935854 K a − ( 3002.5 ⋅1001.1452 ) ) 3002.5
b1 = −
1 (8.2018905K a − 3005938.463) 3002.5
b1 = 0,002731687K a − 1001,1452 22
c1 = −
3002.5 1.0935854 K a 1 0 ( -0,002731687K a + 1001,1452 ) 0,002731687K a − 1001,1452
c1 = −
1 ( 0 − ( 0,002731687K a − 1001,1452 ) ⋅1.0935854 K a ) ( -0,002731687K a + 1001,1452 )
c1 =
( -0,00298733K a + 1094.837774 ) K a ( -0,002731687K a +1001,1452 ) s3 s2
7.5 3002.5
-0,002731687K a +1001,1452 s1 ( -0,00298733K a + 1094.837774 ) K a s0 ( -0,002731687K a +1001,1452 )
1001.1452 1.0935854 K a 0 0
Como se comentó anteriormente los polos positivos coinciden con el número de cambios de signo en la primera columna de la matriz, por lo que ahora ahora se puede ver, que K la estabilidad del sistema depende del parámetro a . D De modo que para determinar la estabilidad del sistema hemos de encontrar para qué rango de valores de K a para los
valores b1 y c1 hacen el sistema estable, resolviendo las ecuaciones obtenidas y encontrando el rango que hace que ambos coeficientes de la matriz sean positivos. Primer parámetro:
( -0,00298733Ka + 1094.837774) Ka ( -0,002731687Ka +1001,1452)
= 0 → Ka = 0
-0,00298733K a 2 1094.837774 K a + =0 +1001,1452 ) ( -0,002731687Ka +1001,1452) ( -0,002731687Ka +1001,1452 − ( -0,002731687K a +1001,1452 ) 0,00298733K a 2 = − ( -0,002731687 -0,002731687K a +1001,1452 )1094.837774 1094.837774K Ka 0,00298733K a 2 0,00298733 = 1 → Ka = 1 → K a = 366493.749937 1094.837774 1094.837774 Ka Por lo que según este factor K a ha de estar entre 0 y 366493.749937 366493.749937. 49937
23
Veamos a continuación que pasa con el otro parámetro: Veamos
-0.002731687K a +1001.1452=0 -0.002731687K a =-1001.1452 Ka =
1001.1452 = 366493.379366 0.002731687
Este segundo valor es más restrictivo que el anterior por lo que habremos de quedarnos con él. El significado de estos valores es que si K a fuese superior a uno de ellos ellos,, el término que corresponde con esta ecuación en la tabla de R Roouth uth- Hurwirtz sería negativo y por lo tanto tendríamos tendríamos un cambio de signo, que implica un polo positivo, haciendo inestable el sistema. Por lo tato el rango de este término será [−∞,366493.379366] , como se ha de cumplir que ambos sean negativos, tenemos que quedarnos con la intersección de ambos intervalos, por lo tanto el intervalo final en el que hemos de trabajar es:
[0,366493.749937] ∩ [−∞,366493.379366] = [0,3664 [0,366493.379366] 93.379366] Mediante la herramienta CC, podemos comprobar esto de una manera rápida y sencilla, con el comando CC>>routh,g. Según el criterio de Routh en esta her herramienta ramienta el sistema es estable dependiendo del parámetro K a en el intervalo [ 0,366493.4] por lo que podemos comprobar que nuestros cálculos anteriores eran correctos. Ahora vamos a analizar si podemos cumplir todos lo loss requisitos del sistema, tan solo con el controlador proporcional que tenemos hasta el momento. Para ello tenemos que recordar las características demandadas a nuestro sistema: ev < 0.1% tr < 0.1s Mp ≤ 30%
Como observamos en los requisitos del sistema, se se nos pide que el error en velocidad no sea superior a 0.1 %, lo que significa que bajo una entrada de escalón unitario ev < 0.001 (0.1% de 1 =0.001). Como hemos obtenido antes, tenemos que este error depende de K a por lo que podemos comprobar si el valor necesario de K a para cumplir este requisito sigue manteniendo el sistema estable:
ev =
915,470525 915,470525 → Ka = Ka ev
Ka =
24
915,470525 = 915470.525 0.001
Por lo tanto, no podemos ajustar el valor de K a para qque ue se cumpla este requisito ya que el límite de estabilidad para un controlador proporcional en nuestro caso está en 366493.379366 según el valor que hemos hallado con el criterio de Routh. Por lo tanto no es suficiente con el controlador proporcional proporcional y tendremos que ir a buscar soluciones más elaboradas.
25
4. Diseño de sistemas P PD D A continuación, vamos a ver cómo realizar el diseño de un controlador PD paara nuestro sistema, y de esta forma ver si somos capaces de estabilizar el sistema con él y a la vez cumplir con los requisitos que necesitamos.
Gb ( s) =
1.0935854 K a s (1001.1452 + 3002.5s + 7.5s 2 )
Los métodos que vamos a aplicar para el sistema, se basan en los sistemas de segundo orden. El orden de un sistema es el grado del exponente de mayor grado de su denominador. N denominador. Nuestro uestro sistema es un sistema de tercer orden orden,, por este motivo nos vemos obligados a obtener un sistema equivalente a este, pero de segundo orden para poder realizar el estudio deseado. Para esto, usaremos la función de transferencia con su denominador descompuesto en producto de factores:
Gb ( s ) =
0.145811K a s ( s + 0.337154 )( s + 399.999618 )
Aquí vemos que tenemos un sistema con un polo muy alejado del centro de coordenadas y del resto de polos. Esto significa que tiene una pequeña aportación y por lo tanto seguramente podamos obviarlo y analizar el sistema como si fuese de segundo orden orden.. Por lo tanto vamos a ver qué ocurre si suponemos:
Gb ( s ) =
0.145811K a 0.145811K a ≈ s ( s + 0.337154 )(( s + 399.999618 ) s ( s + 0.337154 )
Vamos a comparar las respuestas de ambos sistemas en CC, para ver qué ocurre. Introduci Introducimos mos ambas funciones de transferencia, suponiendo Ka = 1 en CC, como G y Geq, ejecutamos las siguientes órdenes para poder comparar ambas graficas: CC>>time,Geq CC>>time,G eq,1,auto ,1,auto a 1,1000 a 1 G [enter] [enter]
26
//Muestra la simulaci simulación ón de G Geq.. //Ajusta el eje temporal para poder compara ambas funciones. //Incluye G en la grafica //Incluye grafica.
En azul podemos ver la señal Geq ( s ) y en rojo Gb (s) . Como podemos observar, los sistemas se diferencian mucho. Esto se debe a que el polo que hemos eliminad eliminado, o, tiene una alta contribución de ganancia al sistema, y por lo tanto no podemos eliminarlo de esta manera. manera. Sin embargo si pasamos esa parte de la contribución del polo a la ganancia del sistema, sistema, la influencia del polo disminuye enormemente de la siguiente forma:
0.145811K a 0.000364528K a 399.999618 G1 ( s ) = = 1 ( s + 399.999618) s ( s + 0.337154 ) s + 1 s ( s + 0.337154 ) 399.999618 399.999618 0.000364528K a 0.000364528K a G1 ( s ) = ≈ s ( s + 0.337154 )( 0.0025s + 1) s ( s + 0.337154 ) Con esto esto hemos pasado a la ganancia del sistema, la contribución del término independiente del polo. Ahora el polo si puede ser eliminado para la realización del análisis, ya que el sistema resultante es muy parecido al anterior. Sin embargo el sistema real, seguirá seguirá siendo un sistema de tercer orden y como se verá más adelante eso nos dará como resultado una mejor respuesta que si fuese de primer o segundo orden. Por este motivo no nos interesa anular este polo definitivamente, definitivamente, aunque el análisis de sistemas de segundo orden es mucho más sencillo. sencillo. Podemos ver a continuación la grafica que representa ambas funciones y en la que se puede apreciar que son prácticamente equivalentes.
27
Un controlador PD consiste en un controlado controladorr con una parte proporcional, al que se le suma un aporte diferencial K p + K D s , y con esto se pretende realizar el control del sistema, mejorando el estado estacionario. Dado un sistema como el siguiente:
θ R ( s)
+
K p + KDs
Va (s)
− θ (s)
ωn 2 s( s + 2δωn )
θ S ( s)
S
La forma habitual de la función de transferencia de un sistema de segundo orden en
ωn 2 . Si analizamos el sistema anterior en bucle cerrado s 2 + 2δωn s + ωn 2 podemos calcular una δ ' que nos de un sistema de segund segundo o orden equivalente pero con bucle cerrado es
un cero más en su función de transferencia, esta equivalencia podemos obtenerla con la siguiente demostración demostración:
Gba ( s )
1 + K D s ) ωn 2 ( = →G s ( s + 2δωn )
bc ( s ) =
(1 + K D s ) ωn 2
(1 + K D s ) ωn 2
s ( s + 2δωn ) s ( s + 2δωn ) = (1 + K D s ) ωn 2 s( s + 2δωn ) + (1 + K D s ) ωn 2 1+ s ( s + 2δωn ) s ( s + 2δωn )
(1 + K D s ) ωn 2 s 2 + 2δωn s + (1 + K D s ) ωn 2
=
(1 + K D s ) ωn 2
s 2 + ( 2δωn + K Dωn 2 ) s + ωn 2 debe equivaler a 2δ 'ωn
2δ ' ωn ≡ ωn ( 2δ + K Dωn ) → 2δ ' ωn ≡ ωn ( 2δ + K Dωn ) → δ ' =
2δ + K Dωn 2
Añadir el controlador PD por tanto, tanto, equivale a añadir un cero al sistema y modificar el coeficiente de amortiguamiento del sistema sistema.. Lo importante es añadir este cero correctamente para conseguir los resultados que deseamos. Es como tener un sistema prototipo de segundo orden donde δ lo hemos aumentado, lo que nos da un sistema menos oscilante. De aquí y dado nuestro sistema equivalente de segundo orden tenemos las siguientes ecuaciones:
Geq ( s) =
0.000364528K a s ( s + 0.337154 )
ωn 2 = 0.000364528K a 2ωnδ = 0.337154
ωn Ahora tenemos 3 incógnitas por determinar: K a δ
28
Determinamos estas estas incógnitas, según el valor que necesitemos dar a δ ' con el fin de que nuestro sistema cumpla las especificaciones especificaciones.. A saber:
Mp =e
δ 'π
−
1−δ '2
≤ 30% tr =
π − arccos δ ' ωn 1 − δ '2
≤ 0.1s
Por tanto hemos de deducir el valor de δ ' que nos resuelva el problema: Mp =e
−
δ 'π 1−δ '2
− ≤ 30% → ln e
δ 'π 1−δ '2
δ 'π ≤ ln ( 0.3) → − ≤ ln ( 0.3) 2 1 − δ '
−δ ' π ≤ -1,203973 1 − δ '2 → δ ' π ≥ 1,203973 1 − δ '2 → (δ ' 2,609355 ) ≥ 2
(δ ' 2,609355 )
2
≥ 1 − δ '2 → δ '2 7,808734 ≥ 1 → δ ' ≥
(
1 − δ '2
)
2
1 7,808734
δ ' ≥ 0,357857
Al introducir un cero por el controlador PD, habrá cambios dinámicos en el sistema, así que mejor aproximar por encima δ ' = 0.5 . tr =
π − arccos δ ' ωn 1 − δ '
2
π − arccos(0.5) 0.1 1 − δ '2
≤ 0.1s → π − arccos(0.5) ≤ 0.1ωn 1 − δ '2
≤ ωn → ωn ≥
π − arccos(0.5) 0.1 1 − 0.25
=
2.0943951 0.0866025
ωn > 24.184003 Por lo tenemos que actualizar el valor de δ ' para que se cumplan nuestros requisitos requisitos.. ω = 30 Para calcular su valor usaremos n y solo nos falta conocer el δ de nuestro sistema sistema: 900 ωn 2 = 0.000364528 K a → 900 = 0.000364528 K a → = Ka 0.000364528 K a = 2468946.144055 2ωnδ = 0.337154 → δ =
0.337154 0.337154 = 2ωn 2 ⋅ 30
δ = 0.00561923
29
δ '=δ +
ωn K D 2
= 0.00561923 +
30 K D = 0.00561923 + 15 K D 2
0.5 = 0.00561923 + 15 K D → 0.494381=15K D → K D =
0.494381 15
K D = 0.0329587 Vamos a redondear el valor de Ka = 2500000 por lo que nnuestro uestro sistema queda de la siguiente forma:
θ R (s)
+
1 + 0.0329587s 0.0329587
− θ ( s)
Va (s)
0.000364528 2500000 s ( s + 0.337154 )
θ S (s)
S
Ahora nos queda hacer las pruebas pruebas del sistema, para comprobarlo, hay que tener en cuenta que el sistema real tiene un polo más que no debemos olvidar olvidar,, y ese es el sistema que realmente debemos probar. La respuesta del sistema queda de la siguiente forma, como se puede ver una sobre oscilación, de 29. 29.3 3% % por lo que en princi principio pio parece que se cumple nuestro requisito. Y tenemos que el tiempo de subida es igual a −3 −1 tr = t (0.9) − t (0.1) = 0.03521 − 3.521⋅10 = 0.031689s < 0.1s siendo tr ( x) = f (t ) es decir la inversa de la función f (t ) .
Observando estos datos podríamos concluir erróneamente, que hemos conseguido satisfacer los requisitos pero aún hemos de comprobar el funcionamiento del sistema 30
real, ya que aquí estamos obviando uno de los polos. El sistema que realmente tenemos es el sig siguiente, uiente, así que vamos a ver cuál es su respuesta real: θ R (s)
+
−θ
1 + 0.0329587s 0.0329587 S
Va (s)
0.000364528 ⋅ 2500000 s ( s + 0.337154 )((1 + 0.0025s )
θ S ( s)
( s)
A continuación podemos ver la respuesta del sistema ante un escalón unitario:
Como se puede observar en la imagen anterior anterior,, tenemos una sobre oscilación del 33. 33.1% % por lo que no cumple los requisito requisito.. Tenemos el siguiente valor para el tiempo de subida tr = t (0.9) − t (0.1) = 0.03521 − 4.93 ⋅10−3 = 0.03028s < 0.1s por lo que este valor si se cumple de sobra. Es interesante observar que el hecho de añadir el controlador PD, nos permite algo que antes no podíamos conseguir:
31
Como puede verse verse, en los datos obtenidos por el método de Routh Routh,, para el sistema original y el sistema al que se le ha añadido el controlador PD, tenemos que antes existía un límite para Ka mientras que una vvez ez introducido el controlador PD ya no tenemos este límite. límite. Por tanto ahora no existe problema para obtener un valor de Ka que nos ayude a eliminar el error. No obstante tenemos el problema de la sobre oscilación del sistema. ¿¿Qué Qué podemos odemos hacer para solucionarlo? solucionarlo?... Pues podemos aumentar el valor de δ ' por ejemplo δ ' = 0.7 para obtener una menor sobre oscilación. Así sí que tenemos que analizar de nuevo el sistema, sistema, de manera que comenzamos un nuevo ciclo de diseño del sistema:
δ '=δ +
ωn K D 2
= 0.00561923 +
30 K D = 0.00561923 + 15 K D 2
0.7 = 0.00561923 + 15K D → 0.694381=15K D → K D =
0.694381 15
K D = 0.046292 Por lo que el sistema final sería: θ R ( s)
+
1 + 0.046292 0.046292s
− θ (s)
Va (s)
0.000364528 ⋅ 2500000 s ( s + 0.337154 )( ) 1 + 0.0025s )
θ S ( s)
S
Comprobaremos a continuación el resultado de la respuesta del sistema ante una entrada escalón unitario:
32
Podemos ver que el sistema, ahora responde con una sobre oscilación muy aceptable, 23.2 23.2% % y podemos quedarnos con este regulador PD. Además claramente podemos ver en la grafica que la subida completa del sistema se realiza en menos de 0.1 segundos por lo qu quee este requisito también queda perfectamente satisfecho, queda comprobar que el erro errorr en velocidad del sistema está dentro también de los parámetros que se han definido como requisitos del sistema.
ev = lim s s →0
ev = lim s s →0
1 s 2 1 + (1 + 0.046292s)
1 0.000364528 ⋅ 2500000 s ( s + 0.337154)(1 + 0.0025s )
1
1
2
s 1 + (1 + 0.046292s) 1 s→ s 0 s 1+
ev = lim
911.32 s ( s + 0.337154)(1 + 0.0025s )
1 (1 + 0.046292s)911.32 s ( s + 0.337154 )( ) (1 + 0.0025s )
1 1 s →0 s s ( s + 0.337154 ) (1 + 0.0025s ) + (1 + 0.046292s)911.32 s ( s + 0.337154)(1 + 0.0025s )
ev= lim
ev= lim s →0
s ( s + 0.337154)(1 + 0.0025s ) 1 s s ( s + 0.337154)( ) (1 + 0.0025s ) + (1 + 0.046292s)911.32 ev=
0.337154 = 3,69962 ⋅10−4 911.32 ev < 0.001
Podemos concluir que con el controlador PD podríamos cumplir todos los requisitos, de manera que tenemos nuestro primer controlado controladorr funcional funcional:
K a = 2500000 GcPD ( s ) = 2500000 (1 + 0.046292 ) G pd ( s ) = 1 + 0.046292s
θ R ( s)
+
1 + 0.046292 0.046292s
− θ (s) S
33
Va (s)
0.000364528 ⋅ 2500000 s ( s + 0.337154 )( ) 1 + 0.0025s )
θ S ( s)
5. Diseño de sistemas PI A continuación, vamos a ver como diseñar un controlador PI P para controlar nuestro sistema, y de esta forma, ver si somos capaces de estabilizar el sistema y a la vez cumplir con los requisitos que necesitamos.
G0 ( s) =
1.0935854 K a s (1001.1452 + 3002.5s + 7.5s 2 )
Al igual que en el caso del controlador PD los métodos que vamos a aplicar para el sistema, se basan en los sistemas de segundo orden. El orden de un sistema es el grado del exponente de su denominador. Nuestro sistema es un sistema de tercer orden, por este motivo nos vemos obligados a obtener un sistema equivalente a este, pero de segundo orden para poder realizar el estudio deseado. Para esto, usaremos la función de transferencia con su denominador descompuesto en producto de factores:
G1 ( s) =
0.000364528K a s ( s + 0.337154 )( ) 0.0025s + 1)
Un controlador PI P consiste en un controlado controlador con una parte proporcional, al que se le KI suma un aporte integral K p + , y con esto se pretende realizar el control del sistema, s mejorando el estado estacionario. Dado un sistema como el siguiente en el que K p = 1 :
θ R ( s)
+
K 11++ I s
− θ (s)
Va (s)
ωn 2 s( s + 2δωn )
θ S ( s)
S
Añadir el controlador P PI representado, equivale a añadir un polo al sistema. Esto tiene el problema, de que afectará a su estabilidad. Aunque Aunque en nuestro caso al convertir el sistema en un sistema de de tipo 2 nos asegura un error de velocidad nulo ev = 0 . Como Como el controlado P PII no nos permite especificar unos valores que tengan un resultado concreto, hemos de probar para ver en qué valores de K I hacen el sistema se mantiene estable. Suponemos K I = 0.5 → Gb ( s ) =
34
( s + 0.5)0.000364528 K a s ( s + 0.337154 )( 0.0025s + 1) 2
El método de Routh nos dice que para este valor de K I , el sistema es inestable para todo valor de K a . Suponemos K I = 0.01 → Gb ( s ) =
( s + 0.01)0.000364528 K a s ( s + 0.337154 )( 0.0025s + 1) 2
Con este valor de K I tenemos que el sistema es estable para los valores
0 ≥ K I ≥ 359282.5 por lo tanto podemos probar para diferentes valores de K a teniendo en cuenta este valor de K I . M p = 77% Ka = 10000 → tr = 0.80986 − 0.24648 = 0.56338s
M p = 90% Ka = 50000 → tr = 0.35211 − 0.10563 = 0.24648s
35
De estos resultados podemos deducir, que conforme aumentemos el valor de K I conseguiremos rebajar el tiempo de subida, pero sin embargo aumentaremos la sobre oscilación y viceversa viceversa.. Dado que por tanto hemos llegado a dos extremos que no satisfacen ninguno de los dos. El que empeoran ya ha sido rebasado y el que mejoran aún no ha conseguido llegar a valores aceptables, concluimos que con este valor de K I no podemos solucionar el problema. Suponemos K I = 0.005 → Gb ( s ) =
( s + 0.005)0.000364528 K a s ( s + 0.337154 ) ( 0.0025s + 1) 2
Con este valor de K I tenemos que el sistema es estable para los valores
0 ≥ K I ≥ 359282.5 por lo tanto podemos probar para diferentes valores de K a teniendo en cuenta este valor de K I . M p = 76.6% Ka = 10000 → tr = 0.80986 − 0.24648 = 0.56338s
M p = 89.9% Ka = 50000 → tr = 0.36972 − 0.10563 = 0.26409s
36
Basándonos en lo que podemos observar de los resultados, parece ser que al rebajar K I rebajamos el valor de sobre oscilación, pero aumentamos el tiempo de subida. Además en los dos casos que hemos analizado hemos llegado a la conclusión de que no se podía mejorar el sistema de forma sustancial, al menos no como para alcanzar los valores qque ue deseamos. Y se puede deducir que no existe valor alguno de K I para el que podamos hallar un K a que resuelva los requerimientos requerimientos.. Además también podemos observar que los sistemas resultantes de un controlado controladorr de este tipo, tienden a tener un alto grado de oscilación, esto seguramente se deba a una combinación del tipo y del grado del sistema obtenido, obtenido, lo que combinado seguramente con otro tipo de reguladores podría dar resultados muy buenos, pero es complicado complicado por si solos.
37
6. Diseño de un controlador respuesta en frecuencia.
PID
empleando
la
En este apartado vamos a realizar el estudio del sistema en el dominio de la frecuencia, y a partir de este análisis diseñaremos un controlador PID. Para ello utiliz utilizaremos aremos herramientas que nos permitirán el análisis de la respuesta en frecuencia del sistema, como los diagramas de Bode de fase y magnitud. Implementaremos este PID como una combinación de una red de adelanto y otra de atraso. Partimos del sistema que estamos estamos analizando: θ R ( s)
+
Va (s)
− θ (s)
0.000364528 ⋅ K a s ( s + 0.337154 )(1 + 0.0025s )
θ S (s)
S
M p ≤ 30% tr = 0.1s ev = 0.001 Para mayor comodidad pasamos la ganancia de los polos al numerador: 0.000364528K a 0.001081191K a 0.337154 G2 ( s ) = = s ( s + 0.337154 ) s 2.966004s + 1)(( 0.0025s + 1) ( 0.0025s + 1) ( 0.337154 El diseño de un regulador es un proceso de tanteo, debemos ajustar el regulador siguiendo los siguientes pasos de forma iterativa hasta obtener el resultado deseado. Calcular el error en estado estacionario y determinar K a Intentar ajustar un PD Ajustar un PI
38
6.1 Calcular el error de velocidad (determinar K a ) ev = lim sE (s) = lim s s →0
ev = lim s s →0
ev = lim s →0
ev = lim s 0 s→
ev = lim s→0
s →0
R( s ) 1 + Gp ( s)
1 1 0.000364528 ⋅ Ka s 21 11+ + s ( s + 0.337154)(1 + 0.0025s ) 1 1 0.001081191⋅ Ka s 1+ s ( 2.966004s + 1)( 0.0025s + 1)
1 1 s ( 2.966004s + 1)( 0.0025s + 1) + 0.001081191⋅ Ka s s ( 2.966004s + 1)( ) ( 0.0025s + 1)
s ( 2.966004s + 1)( 0.0025s + 1)
1 1 = s ( 2.966004s + 1)( 0.0025s + 1) + 0.001081191⋅ Ka s 0.001081191⋅ Ka ev =
924.905960 Ka
Por lo que el valor de K a necesario para conseguir el error en velocidad que nos piden los requisitos del sistema es: es
ev =
924.90596 924.90596 → ev = ≤ 0.001 Ka Ka 924.90596 ≤ Ka 0.001 Ka ≥ 924905.96
Ahora tenemos que comprobar si el sistema es estable para este valor de K a mediante el criterio de Routh.
39
Según el software CC, el sistema será estable para valores de Ka ∈ [ 0,370274.2] lo que indica que no podemos eliminar el error en velocidad por completo de esta manera manera,, ya que para eso se necesita un valor de K a ≥ 924905.96 . Por tanto sacrificaremos de momento el error y nos centraremos en el análisis de la estabilidad estabilidad.. Analizaremos Analizaremos que ganancia nos interesa aplicar en fu función nción del tiempo de subida que queremos para el sistema. Podemos ver en el diagrama, que el punto de cruce de magnitud, se encuentra en frecuencias muy bajas. Esto significa que el sistema, es un sistema muy lento, por lo que habremos de usar un K a alto para ajusta ajustarr el sistema. Para decidir este valor veremos en qué valores necesitamos el cruce de magnitud para cumplir con el tiempo de subida del sistema.
Para los sistemas de segundo orden tenemos dos aproximaciones que podemos emplear para realizar modificaciones en su comportamiento. En concreto tenemos que para 0 < δ < 0.6 podemos aplicar: aplicar
1 ω PCM ≈ t r 0 < δ < 0.6 δ ≈ MF 100 Como nuestro sistema se comporta casi como si fuese un sistema de segundo orden podemos emplear estas ecuaciones para buscar los parámetros deseados del sistema siempre y cuando δ esté dentro del rango mencionado mencionado.. De este modo, tenemos que hallar el valor de δ para que nuestro sistema cumpla las especificaciones y ver si podemos cumplirlas con un valor dentro del rango. rango Para obtener este valor, usamos como siempre la formula: 40
Mp =e
−
δ 'π 1−δ '2
≤ 30%
Por tanto tenemos:
Mp = e
−
δπ 1−δ 2
− ≤ 30% → ln e
δπ 1−δ 2
δπ ≤ ln ( 0.3) → − ≤ ln ( 0.3) 2 1 − δ
−δπ ≤ -1,203973 1 − δ 2 → δπ ≥ 1,203973 1 − δ 2 → ( 2,609355δ ) ≥ 2
( 2,609355δ )
2
≥ 1 − δ 2 → δ 2 7,808734 ≥ 1 → δ ≥
(
1− δ 2
)
2
1 7,808734
δ ≥ 0,357857
De modo que sí s podemos cumplir con las especificaciones del sistema cumpliendo con la restricción para aplicar las aproximaciones. aproximaciones. Podemos por tanto, obtener la frecuencia a la que queremos transportar nuestro punto de cruce de m magnitud agnitud mediante la siguiente fórmula:
ω PCM =
1 1 → ω PCM = = 10(rad / sec) tr 0.1
Veamos a continuación cual es la ganancia que necesitamos para que nuestro punto de cruce pase pase a la frecuencia de 10 rad/s rad/s.. En el diagrama de Bode podemos ver que la ganancia que necesitar necesitaríamos íamos es de 108 dB aproximadamente.
Sin embargo, la red de adelanto produce una atenuación de la magnitud, de manera que el nuevo punto de cruce de magnitud tras aplicar la red queda centrado en el punto 41
del diagrama de Bode actual cuya magnitud es de −10log α . Por tanto hemos de emplear una ganancia más alta aunque esta haga al sistema inestable. inestable. Más adelante podremos solucionar esto con un adelanto de fase. Por tanto probaremos con 25 25 rad/s. rad/s
Vemos que para llevar nuestro PCM a frecuencias de 25 5 rad/s tenemos que aportar 125 25 dB de ganancia, por tanto tenemos que K a debería valer: 125 125 = 20 log K a → = log x → 106.25 → K a ≥ 1778279.410039 20 Probaremos con una ganancia de 2000000 00000 para curarnos en salud salud.. V Vamos amos a ver a continuación como se ha modificado el diagrama de Bode Bode.
Como podemos ver ver,, tenemos un margen de aproximadamente unos 17 dB para añ añadir adir la red de adelanto de fase. Hay que tener en cuenta que el valor estimad estimado o de la 42
frecuencia de 10 rad/s rad/s,, que nos permite conseguir el tiempo de subida de 0.1s, 0.1s es aproximada. No obstante 17 dB sigue siendo un buen margen. Se puede comprobar que aproximada. el sistema será inestable ya que en el punto de cruce de magnitud, la fase está por debajo de -180 180 180º. Tenemos entonces, la siguiente función de transferencia: G3 ( s ) =
2162.382 s ( 2.966004 s + 1)( ) 0.0025s + 1)
Y por por último nuestro sistema queda de la siguiente forma:
θ R ( s)
+
Va (s)
− θ (s) S
43
2162.382 s ( 2.966004s + 1)(( 0.0025s + 1)
θ S (s)
6.2 Diseño de la red de avance de fase En primer lugar veamos la respuesta de nuestro sistem sistemaa ante una entrada escalón:
Cómo podemos ver el sistema ofrece un unaa respuesta muy oscilante y no parece responde bien en régimen estacionario. A continuación vamos a ver el diagrama de Bode para analizar el comportamiento del sistema sistema.. Para esto introducimos introducimos los comandos: CC>>siggy,gfk,6,a CC>> siggy,gfk,6,a //Podemos Podemos ver el diagrama de Bode del sistema (fase,magnitud) D 3,2 [enter] [enter] //Muestra Muestra la gráfica en dB
Podemos ver que la fase en el punto de cruce de magnitud es de -183.08 183.08 lo que significa que tenemos un sistema inestable como comentamos en el apartado anterior. Esto explica el comportamiento del sistema ante una entrad escalón. El Margen de fase de nuestro sistema es por tanto 0. 0 El margen de fase representa, el valor máximo para el 44
cual se puede asegurar que un cambio de fase en el punto de corte de magnitud no hace inestable el sistema. De este modo si es muy pequeño el sistema se va acercando a un sistema críticamente estable y por este motivo aumenta sus oscilaciones. La solución pasa por aumentar el m margen argen de fase del sistema, ya que esto nos modifica el régimen transitorio. Un sistema es estable siempre que la fase en el punto de cruce de magnitud, esté comprendida entre 0 y -180º -180º.. Por tanto tenemos que hacer que el sistema entre dentro de este rango en el PCM y además con un margen de fase suficientemente alto para que se cumpla nuestro requisito de sobre oscilación oscilación. Para realizar este tipo de modificaciones en el diagrama de Bode, tenemos las redes de adelanto y atraso de fase. Estas redes, nos per permiten miten cambiar la fase en un punto concreto del diagrama consiguiendo un resultado aproximado de lo que queremos obtener obtener. Como hemos visto en el apartado anterior necesitamos cumplir que δ ≥ 0,357857 para que la sobre oscilación no supere el 330%, 0%, así que tomaremos valor, δ = 0.5 ya que estas aproximaciones son más certeras para 0 > δ > 0,6 . Con este valor sobreestimamos un poco, poco para cumplir nuestro objetivo y así nos aseguramos un margen alto de tolerancia tolerancia. Por último para obtener el M F que necesitamos calculamos δ =
MF 100
de modo que tenemos tenemos:
MF → M F = 100δ = 100 ⋅ 0.5 = 50º 100 usaremos 55º por si acaso ya que el sist sistema ema
δ=
no es realmente de 2º orden y sobre osci oscilará lará un poco más de lo previsto
φad = 55º Conocido el φad que quere queremos os aumentar, podemos calcular una función en forma de campana de Gauss, tal que tenga una gran aportación en el punto de corte de magnitud de nuestro sistema y este aporte decrezca en los alrededores alrededores,, intentando minimizar las modificaciones en el resto del diagrama de Bode. Este punto se encuentra en la frecuencia cuya magnitud se corresponde corresponderá rá con −10log α . A su vez esta frecuencia se puede expresar de la forma adelanto en
45
10−1 102 y . T T
1 T α
quedando por tanto los extremos de la red de
La red de adelanto queda del siguiente modo: K → Ganancia adicional del sistema. Permitir Permitiráá ajustar el diagrama de c magnitud para aprovechar el punto justo de máximo aporte de fase de la red. T →T = 1 Factor que nos da la frecuencia de inicio y fin ωm α 1 + Ts Ga ( s ) = K cα 1 + α Ts 10-1 102 de aportación de la red y . T T α → α = 1 − senφad Aporte de ganancia en el punto donde qu quedará edará 1 + senφad el nuevo PCM.
El diagrama de bode de la red de avance de fase quedaría de la siguiente forma:
De aquí determinamos α para el avance de fase que deseamos como:
α=
46
1 − senφad 1 − sen(55º ) 0,180848 = = = 0,099413 1 + senφad 1 + sen(55º ) 1,819152
Sabiendo que ωm se corresponde con la frecuencia en el diagrama de Bode para el valor de magnitud −10log α .
G( jω ) ω =-10logα (dB)=-10log0.099413=10.0255682 (dB)=-10log0.099413=10.0255682(dB) (dB) m
Una vez calculado el valor habrá que buscar en el diagrama de Bode este valor y su correspondiente frecuencia en radianes por segundo. segundo. Este punto será el nuevo PCM de del sistema y por tanto donde nos interesa centrar la red de adelanto, para aprovechar al máximo la aportación de fase fase..
Del diagrama podemos extraer que la frecuencia en la que se realizará el máximo apo aporte rte de fase ha de ser 15 15 rad/s. rad/s. E Esto sto es una buena señal ya que esta frecuencia parece permitir que la respuesta del sistema sea lo suficienteme suficientemente nte rápida como para obtener un tiempo de respuesta aceptable aceptable. Por tanto ωm = 15 rad/seg . T=
1 1 = 15 0.099413 4.729474 T = 0.21144
Finalmente tenemos una red de avance de fase igual a: Ga ( s ) = 0.099415
47
1 + 0.21144s 1 + 0.0210199s
Como podemos observar en el diagrama de Bode, la red de adelanto ha modificado ligeramente la frecuencia del punto de corte de magnitud, que antes estaba en torno a los 26.63 rad/s y ahora se encuentra en 15.1 rad/seg.
Sin embargo, la red de adelanto se ha centrado prácticamente en el punto de corte de magnitud. El nuevo margen de fase del sistema es 180º −125.88º = 54.12º ≃ 55º +0.07º que queríamos conseguir. Podemos ver a continuación la respuesta del sistema ante una entrada escalón, escalón, tenemos tenemos un sistema que cumple los requisitos de sobre oscilación y de tiempo de subida subida.
Como podemos ver el margen que tenemos para la sobre oscilación es grande grande.. Podríamos aplicar una ganancia al sistema sistema, si queremos queremos, para hacerlo un poco más rápido rápido, sin n comprometer demasiado la sobre oscilación que que tenemos en el 22.4%. Sin embargo como la red de retraso no introducirá perjuicio en el sistema y ya se cumplen los requisitos podemos dejarlo así.
48
6.3 Diseño de la rred ed de retraso retraso de fase Llegados a este punto, partiremos de nuestra función completa en este momento, que incluye, la ganancia K p y la red de adelanto de fase. Con ella analizaremos las bajas frecuencias del diagrama de Bode que corresponden con la respue respuesta sta del sistema en estado estacionario y por tanto aquí podremos solucionar nuestro problema con el error de velocidad. G4 ( s ) = 0.099415
1 + 0.21144s 2162.382 1 + 0.0210199s s ( 2.966004 s + 1)(( 0.0025s + 1)
El error en velocidad actual del sistema es igual a:
ess = lim s s →0
ess = lim s s →0
R( s ) 1 + Gb 4 (s) 1
1 + 0.099415
1 + 0.21144s 2162.382 1 + 0.0210199s s ( 2.966004s + 1)( ) ( 0.0025s + 1)
1 s 21
1 1 s→0 214.973207(1 + 0.21144s) s 1+ (1 + 0.0210199s ) s ( 2.966004s + 1)( ) ( 0.0025s + 1)
ev = lim
ess = lim s →0
1 1 .21144s) s (1 + 0.0210199s ) s ( 2.966004s + 1))(( 0.0025s + 1) + 214.973207(1 + 00.21144 (1 + 0.0210199s ) s ( 2.966004s + 1)( 0.0025s + 1)
(1 + 0.0210199s ) s ( 2.966004s + 1)( 0.0025s + 1) 1 s →0 (1 + 0.0210199s ) s ( 2.966004s + 1) .21144s) s )(( 0.0025s + 1) + 214.973207(1 + 00.21144
ess = lim
ess = lim s→0
2.966004s + 1) ( 0.0025s 0.0025 s + 1) (1 + 0.0210199s ))(( 2.966004
0.21144 ) (1 + 0.0210199s ) s ( 2.966004s + 1)( 0.0025s + 1) + 214.973207(1 + 0.21144s ess =
1 214.973207
ess = 0.00465174
49
Con esta red de atraso, podemos rebajar la fase en las bajas frecuencias, consiguiendo mejorar el error en estado estacionario del sistema. Dado que la red de 102 atraso deja de afectar al resto del diagrama en las frecuencias ωm = si la colocamos T en la frecuenc frecuencia ia de cruce, siempre que se cumpla que 2 4 2 ωPCM 10 ω 10 10 < →T < →T < , siendo PCM = ωm tendremos una red de atraso 100 T ωPCM ωm 100 de fase que no afectará a la fase en el punto de cruce de magnitud, por lo que no cambiará la sobre oscilación actual de nuestro sistema.
1 1 → Kv = = 214.973322 , comenzamos la construcción Kv 0.00465174 de la red de retardo: 1 + Ts Gr ( s) = β 1+ T β s
Dado que ess =
De forma que se cumpla que Kv = 214.973322β y como sabemos que queremos 1000 Kv ≥ 1000 tenemos que 214.973322β ≥ 1000 → β ≥ → β ≥ 4.6514 214.973322 Tomamos por ejemplo β = 4.8 . L Lo o tomamos un poco más grande por si acaso.
50
Ahora vamos a ver el diagrama de Bode de nuestra función de transferencia y en ella buscamos la frecuencia de cruce de magnitud.
Dichaa frecuencia se halla en 15.1 rad/s la tendremos en cuenta para intentar no Dich estropear, o estropear lo menos posible la parte del sistema que ya hemos diseñado diseñado.. Esto lo conseguiremos alejando la red de retraso, teniendo en cuenta la ecuación comentada anteriormente, si usamos una frecuencia frecuencia de máxima perdida de fase ω 15.1 ωm = PCM = = 0.151 rad/seg tendremos una red de adelanto de fase que no nos 100 100 100 perjudicará siempre y cuando se cumpla que T < . Por tanto buscamos T: ωm
T=
1
ωm β
→T =
1 1 = 0.151 4.8 0.330824
T = 3.0227553 <<<
100 = 662.25 0.151
De forma que la red de retardo de fase queda como sigue sigue,, y no afectará al resto del diseño diseño: Gr ( s ) = 4.8
51
1 + 3.0227553s 1 + 14.509225s 14.509225
Tras aplicar la red de retraso construida, podemos ver que el margen de fase de nuestro sistema no se ha modificado prácticamente prácticamente,, ahora tenemos 180º −126.88º = 53.12º ≈ 54.12º por lo que lo más probable es que hayamos aumentado un poco la sobre oscilación del sistema.
Se ha rebajado el tiempo de subida a 0. 0.955 955 , aunque ha aumentado un poco la sobre 955s, oscilación hasta 24.04 24.04% como nos imaginábamos imaginábamos.
La función de transferencia del sistema completo queda como: G5 ( s ) = 4.8
1 + 3.0227553s 1 + 0.21144s 2162.382 0.099415 1 + 14.509225s 1 + 0.0210199s s ( 2.966004 s + 1) ( 0.0025s + 1)
Compararemos a continuación el diagrama de Bode anterior al nuevo, para ver que en las altas frecuencias se ha mantenido la respuesta en frecuencia del sistema. En azul y Rojo podemos ver la magnitud y fase del nuevo sistema respectivamente, respectivamente, mientras que en verde y negro podemos ver el sistema antiguo. Podemos fijarnos que en las zonas cercanas al punto de cruce de magnitud, magnitud el diagrama se hhaa mantenido prácticamente igual. All contrario pasa en las bajas frecuencias. Vemos que la red de prácticamente atraso de fase ha rebajado la fase en torno a la frecuencia 0.15 rad/s rad/s, como habíamos 52
diseñado y que se produce un pequeño aporte de ganancia en las bajas frecuencias del sistema lo que mejora el régimen estacionario y por tanto ayudará a reducir su error.
A continuación comprobaremos que el error en estado estacionario para la velocidad del sistema se cumple: G5 ( s ) = 4.8
1 + 3.0227553s 1 + 0.21144s 2162.382 0.099415 1 + 14.509225s 1 + 0.0210199s s ( 2.966004s + 1)(( 0.0025s + 1)
K v = lim sG ( s ) = s 4.8 s →0
K v = lim s 4.8 ss→ →0
1 + 3.0227553s 1 + 0.21144s 2162.382 0.099415 1 + 14.509225s 1 + 0.0210199s s ( 2.966004s + 1)( 0.0025s + 1)
1 + 3.0227553s 1 + 0.21144s 2162.382 0.099415 1 + 14.509225s 1 + 0.0210199s s ( 2.966004s + 1)( 0.0025s + 1)
(1 + 3.0227553s )(1 + 0.21144s )1031.871391 s→0 (1 + 14.509225s ) )((1 + 0.0210199s )( ) ( 2.966004s + 1)( 0.0025s + 1)
K v = lim
K v = 1031.871391 ev =
53
1 1 = = 0.000969113 K v 1031.871391
Por último analizaremos el sistema resu resultante ltante en bucle cerrado cerrado,, con el comando “stability” de CC. De esta manera podemos obtener los polos de la función en bucle cerrado cerrado, y observar en la imagen a continuación los polos de la función así como la función de transferencia en bucle cerrado.
Como análisis análisis, podemos empezar comentando que ten tenemos emos 2 polos conjugados puros que serán el eje principal del funcionamiento del sistema. U Un n polo en -402 402 que no tendrá a penas peso en la función ya que se encuentra muy alejado de los polos conjugados puros puros.. Por otro lado nos quedan otros dos polos en -0. 0.3308 33083236... 33083236... y otro en -88.435 435... ..... Como comprobaremos a continuación, estos polos deber deberían ían cancelarse por los ceros, al menos en gran medida. G5 ( s ) = 4.8
1 + 3.0227553s 1 + 0.21144s 2162.382 0.099415 1 + 14.509225s 1 + 0.0210199s s ( 2.966004 s + 1) ( 0.0025s + 1)
1 = −0.330824 ≈ −0.33083236... 3.0227553 1 1 + 0.21144s =0 → s== −4.730369 ≈ −8.43507... 0.21144 1 + 3.0227553s=0 → s=-
Podemos ve verr que el polo -0.3308 - 3308 es cancelado por el cero -0. 0.330824 330824, 330824, mientras que parte del polo en -8.43507 8.43507 es cancelado por el cero en -4.730369 4.730369 4.730369.. Sin embargo este último no es cancelado del todo, lo que deja una pequeña parte del polo contribuyendo al sistema. Esto consigue que el sistema se asemeje mucho al comportamiento de un sistema de segundo orden pero sin llegar a ser un sistema de segundo orden puro. No obstante el mismo hecho de que sea de tercer orden, permite que el sistema sea menos oscilante a mayores vvelocidades elocidades elocidades,, y por tanto que podamos obtener tan buenos resultados con él. Hemos obtenido por tanto nuestro segundo regulador satisfactorio:
54
K a = 2000000 1 + 3.0227553s 1 + 0.21144s 1 + 0.21144s GcPID ( s ) = 4.8 0.099415 Gav ( s ) = 0.099415 1 + 14.509225s 1 + 0.0210199s 1 + 0.0210199s 1 + 3.0227553s Gret ( s ) = 4.8 1 + 14.509225s
55
7. Diseño de un PID por las reglas de Zigher Zigher-Nicols Nicols En un sistema real, en la mayoría de las ocasiones es muy difícil o imposible determinar con exactitud la función de transferencia del sistema que se desea controlar. En estos casos, y con el requisito de no necesitar una respuesta con unas restricciones determinadas (tiempo de respuesta, sobre sobre oscilación máxima, etc.), pueden emplearse las reglas de ZigherZigher Zigher-Nicols. Nicols. Estas reglas nos permiten garantizar la estabilidad del sistema, y el control del mismo mediante un PID General, pero sin poder especificar los parámetros de respuesta del mismo de forma exacta. exacta. A continuación vamos a crear un controlador PID por este método, así podremos ver que el resultado obtenido es mucho pero que el obtenido de forma analítica. analítica θ R ( s)
+
K p (1 +
− θ (s)
1 + K D s) KI s
Va (s)
0.000364528 θ S (s) s ( s + 0.337154 )((1 + 0.0025s )
S
Dado un sistema como el anterior, necesitamos encontr encontrar ar el valor para las constantes K p , K I y K D . Para ello usaremos como hemos comentado las reglas de Zigher Zigher-Nicols. Nicols. En primero lugar necesitamos determinar qué valor necesitamos para K p que haga oscilante el sistema. Para hallar este valor, debemos eliminar el efecto del controlador proporcional y diferencial. En un sistema real, tendríamos que dar valor nulo a la constante K I e infinito a la constante K D . De esta manera podemos modificar el valor de K p buscando uno que haga el sistema críticamente estable, y con él determinar el resto de parámetros. En nuestro caso como conocemos la función de transferencia del sistema, por lo que podemos emplear el método de Routh para determinar este valor limite sin tener que hacer pruebas: pruebas
Encontrar ncontrar este valor es muy sencillo, sin embargo en un caso real, este sería un trabajo de monitorización monitorización y tanteo del ajuste, hasta que se encontrase este valor. Como podemos ver, tenemos el valor 370274.1.
56
A continuación se muestra la respuesta del sistema en el tiempo, con el valor límite para K p . También podemos ver que el periodo periodo de la señal es T=0.54s T=0.54s. 0.54s
Según las reglas de Zigher Zigher-Nicols, Nicols, si nombramos al valor límite para K p cómo K cr y al periodo de la señal en este instante como Pcr , tenemos que el valor que hemos de establecer para los parámetros del PID son: K p = 0.6 ⋅ K cr = 0.6 ⋅ 370274.1 = 22244,46 PID Zigher-Nicols K I = 0.5 ⋅ Pcr = 0.5 ⋅ 0.54 = 0.27 K = 0.125 ⋅ P = 0.125 ⋅ 0.54 = 0.0675 cr D
Podemos observar, que el sistema obtenido es mucho más oscilante que el que obtuvimos mediante las redes de adelanto y atraso. Su sobre oscilación es mayor y es más lento. Sin embargo el sistema es estable y se puede configurar de una forma muy sencilla con sistemas reales. Además puede usar usarse se como base para realizar ajustes más precisos mediante tanteo. tanteo 57
8. Construcción de un PID discreto. En este apartado, vamos a muestrear el sistema, buscando el sistema digital equivalente que represente nuestro proceso, de manera que lo podamos modelar mediante la transformada Z Z,, y a continuación realizar un análisis del mismo, diseñando un PID utilizando el método del lugar de llas as raíces raíces.. Para poder muestrear nuestro sistema, tenemos que emplear emplear un periodo de muestreo que dependerá de la rapidez de nuestro sistema. Buscaremos un periodo de muestreo, que nos permita reconstruir nuestro sistema. Este periodo de muestreo lo estableceremos estableceremos en función del tiempo de subida, determinando una cantidad de muestras para la subida subida. G2 ( s ) =
0.001081191K a s ( 2.966004 s + 1)( ) 0.0025s + 1)
Para poder diseñar el controlador digital, hemos de modelar el sistema como un sistema discreto, más un bloqueador de orden cero cero. Por tanto nuestro sistema y regulador seguirán el siguiente esquema:
θ R (Z ) E(Z )
+
−θ (Z ) S
U (Z )
Gc ( z )
0.001081191K a 1 − eTs U (s) s ( 2.966004s + 1)(( 0.0025s + 1) s
θ S ( s)
Bloqueador de orden cero A/D D
Donde U ( s) equivale a VR (s) del sistema anterior, de modo que la señal de control será la entrada de la planta. Con uun n valor de Ka = 1 podemos ver en el diagrama siguiente siguiente, que el sistema es muy lento, por este motivo será recomendable partir de una ganancia más elevada y así obtener un sistema más rápido de entrada.
58
Usando el criterio de Routh podemos ver que el sistema es estable entre 0 y 370270.2, vamos a coger un valor aleatorio por ejemplo un valor intermedio como 250000. Por tanto Ka = 250000 .
De esta manera nuestra función de transferencia para este valor de ganancia ganancia queda como: G7 ( s ) =
270.2978 s ( 2.966004 s + 1)( ) 0.0025s + 1)
Veamos a continuación su respuesta en el tiempo:
Tenemos un sistema más rápido que el anterior, y ahora hemos de muestrearlo. Para seleccionar el periodo de muestreo Ts tendremos tendremos en cuenta, el tiempo de subida deseado y el de la respuesta actual del sistema sistema.. Como podemos ver en la imagen de la derecha, el sistema tarda en alcanzar el punto de máxima oscilación aproximadamente 0.32s 0.32s, podríamos determinar el tiempo de muestreo en ffunción unción de este valor, por ejemplo tomando tomando 50 muestras como mínimo en este tiempo, tiempo, interpretaremos a esto como una subida en el siguiente análisis: análisis
ωs = ωs =
tr (segundos/subida) 50(muestras/subida)
0.32(segundos/ subida ) 50 = (muestras/segundo)=156 ≈ 200 0.32 50(m/ subid a)
ωs =200 muestras/segundo → Ts =
59
1 = 0.005s 0.005 200
Si esto no fuese suficiente, tendríamos que volver a tomar muestras para calcular de nuevo el sistema. Por Por tanto discretizamos nuestro sistema con este periodo de muestreo muestreo:
Podemos ver a continuación la función de transferencia de nuestra planta discretizada:
4.922663 ⋅10−4 ( z 2 + 2.624847 z + 0.373609) G6 ( Z ) = ( z − 1)( z − 0.996634)( z − 1.831564 ⋅10−2 ) Veamos la respuesta del sistema frente a una entrada escalón:
El resultado oobtenido btenido es un sistema inestable porque no hemos podido reconstruir la señal original. original. E Esto sto es signo de que no hemos cogido las muestras necesarias. Por tanto hemos de usar un periodo de muestreo menor. Como se ha comentado anteriormente el sistema tarda en alcanzar su punto de máxima oscilación 0.32s, sin embargo el sistema que queremos conseguir tiene que tener un tiempo de subida de 0.1s, de modo que usaremos este último para determinar el Ts del sistema ya que nos dará una frecuencia de muestreo mayor. mayor Probaremos esta vez con 100 muestras para la subida de la señal deseada, de modo que tendríamos: 60
ωs = ωs =
tr (segundos/subida) 100(muestras/subida)
0.1(segundos/ subida ) 100 = (muestras/segundo)=1000 0.1 100(m/ subid a)
ωs =1000 muestras/segundo → Ts =
1 = 0.001 0.001s 1000
Si esto no fuese suficiente, de igual modo habría que repetir la operación operación:
Podemos ver a continuación la función de transferencia de nuestra planta discretizada:
G6 ( Z ) =
−1.062405 ⋅10−13 ( z 3 − 5.189153·107 z 2 − 1.883729·108 z − 4.248368·10 4.248368·107 ) ( z − 1)( z − 0.996629)( z − 0.67032)
Veamos la respuesta del sistema frente a una entrada escalón:
61
Esta vez parece que si hemos podido reconstruir la señal correcta, podemos ver que el tiempo en llegar al punto de máxima oscilación es aproximadamente el mismo, 0.33s. Obtengamos ahora su ecuación característica característica: 1 + G ( Z ) = 1 + G6 ( Z ) = 1 +
−1.062405 ⋅10 −13 ( z 3 − 5.189153·107 z 2 − 1.883729·108 z − 4.24836 4.248368·10 8·107 ) ( z − 1)( z − 0.996629)( z − 0.67032)
Buscamos los polos: 1 + G6 ( Z ) = 1 + 1 + G( Z ) =
−1.062405 ⋅10−13 z 3 − 5.51892·10−6 z 2 − 2.00127·10−5 z − 4.51349·10−6 ( z − 1)( z − 0.9996629)( z − 0.67032)
z 2 − 2.669983z 2 + 2.340077 z − 0.6700941 − 1.062405 ⋅10−13 z 3 − 5.518 5.51892·10 92·10−6 z 2 − 2.00127·10−5 z − 4.51349·10−6 ( z − 1)( z − 0.9996629)( z − 0.67032) 1 + G(Z ) =
z 2 − 2.669977 z 2 + 2.340097 z − 0.6700895 =0 ( z − 1)( z − 0.9996629)( z − 0.67032)
Una cosa importante es comentar que el tipo de un sistema en Z equivale a la cantidad de polos en z=1, por lo que comprobamos que nuestro sistema resultante sigue siendo de tipo 1 y tercer orden. Para realizar el diseño de un controlado mediante el lugar de las raíces, primero primero hemos de encontrar los valores de δ y de ωn para que el sistema cumpla las características requeridas. Las siguientes ecuaciones son para sistemas de segundo orden, por lo que esto será un valor aproximado aproximado..
Mp =e
62
−
δ 'π 1−δ 2
≤ 30% tr =
π − arccos δ ωn 1 − δ 2
≤ 0.1s
Mp =e
−
δπ 1−δ 2
− ≤ 30% → ln e
δπ 1−δ 2
δπ ≤ ln ( 0.3) → − ≤ ln ( 0.3) 2 1 − δ
−δπ ≤ -1,203973 1 − δ 2 → δπ ≥ 1,203973 1 − δ 2 → ( 2,609355δ ) ≥ 2
( 2,609355δ )
2
≥ 1 − δ 2 → δ 2 7,808734 ≥ 1 → δ ≥
(
1− δ 2
)
2
1 7,808734
δ ≥ 0,357857 Por ejemplo podemos buscar un valor superior a 0.4 dado que el sistema real no es de segundo orden. tr =
π − arccos δ ωn 1 − δ
2
π − arccos(0.5) 0.1 1 − δ
2
≤ 0.1s → π − arccos(0.5) ≤ 0.1ωn 1 − δ 2
≤ ωn → ωn ≥
π − arccos(0.5) 0.1 1 − 0.25
=
2.0943951 0.0866025
ωn > 24.184003 En este caso de igual modo podemos buscar un valor de ωn entre 25 y 30 para no quedarnos cortos.
Mediante los valores de ωn y δ podemos obtener el valor aproximado de los polos, que harían que el sistema cumpla con nuestras características. Estos polos sin embargo serán los polos conjugados puros para un sistema de segundo orden. Puesto que nuestro sistema es un sistema de tercer orden, tendremos que aproximarnos a ellos y sobrepasar un poco para ajustar el sistema: sistema:
63
K ⋅ ωn 2 K ⋅ ωn 2 G ( s) = 2 = s + 2·δωn ·s + ωn 2 ( s + σ − ωd j )( s + σ + ωd j ) s 2 + 2·δωn ·s + ωn 2 = 0 → s1,2 = −δωn ± ωn 1 − δ 2 j
σ = δ ⋅ ωn
s1,2 = −σ ± ωd j ωd = ωn · 1 − δ 2
Por lo que tenemos:
σ = 0.4 ⋅ 25 ≥ 10 ωd = 0.25· 1 − ( 0.4 )
2
s1,2 ≈ −0.229129 ± 10 j ≥ 0.229129
Sabemos que existe una relación entre los espac espacios ios de convergencia s y Z. Z M Mientras ientras que en el espacio s la región de convergenc convergencia ia de los sistemas se limita a la parte real negativo del espacio, espacio, en Z se transporta a la circunferencia unidad del espacio.
Tenemos pues que la relación existente entre los polos en los diferentes planos es:
z1,2 = eT ·s = eT ( −σ ±ωd j ) = e −Tσ e ±T ωd j z1 = e −T ·σ
−0.01·10 ±0.01·0.221929 j e = e −0.1e ±0.00221929 j z1,2 = e ∠z1 = θ = T ·ωd z1,2 = 0.90484·e ±0.00221929 j Este valor se corresponde en el espacio con los polos conjudaso puros:
z1,2 = 0.90484·e ±0.00221929 j → z1,2 = 0.904838 ± 2.008100715·10−3 j 64
Tenemos por tanto el siguiente lugar de las raíces: raíces
Como podemos extraer de dell diagrama anterior, anterior, nuestros polos están muy lejos de los polos que queremos obtener para cumplir los requisitos del sistema. Por tanto, tenemos que modificar el lugar de las raíces para hacer que este contenga los polos que necesitamos. Par Paraa esto debemos emplear un controlado PD, de manera que con un cero cancelemos uno de los polos, y con el polo añadido modifiquemos el lugar de las raíces para que pase por los polos de nuestros requerimientos. GPD ( z ) = K d
z + zcd z + z pd
En primer lugar hemos de ver el valor del polo que vamos a cancelar:
El polo que nos interesa cancelar es el que se encuentra más cercano al eje de coordenadas ya que si eliminamos el polo en z=1 bajaríamos el tipo del sistema. Podemos odemos ver que en las graficas ddee Mathlab se aproximan los valores, por lo tanto buscaremos el cero en z = 0.9996629 . GPD ( z ) = K d
z − 0.9996629 z + z pd
Ahora tenemos que determinar donde es necesario añadir el nuevo polo, para que el nuevo lugar de las raíces contenga los polo polos que ue hemos definido anteriormente, para ello hemos de utilizar el criterio del argumento argumento. 65
El criterio del argumento establece que para que un punto z pertenezca al lugar de las raíces, debe cumplirse la siguiente relación de argument argumentos: os:
∑ ∠ ( z − c ) − ∑ ∠ ( z − p ) = ( 2r + 1) π i
i
r ∈ℤ
En nuestro caso, tenemos que el punto z en cuestión es el punto ze y que los ceros y los polos son, los ceros y los polos del sistema en bucle abierto sin contar el cero añadido y el polo cancelado ya que se anulan entre sí. Por último se añadirá el que queremos obtener en función de todo los demás demás.
α = ∠ ( ze + zc , p ) = arctan
Im( ze ) ⇒ ze = 0.904838 ± 2.008100715·10 −3 j Re( ze ) − zc , p
CEROS Z c1 = −0.241611
α1 = ∠ ( ze + zc1 ) = arctan
Im( ze ) 2.008100715·10−3 = arctan Re( ze ) − zc1 0.904838 + 0.241611
α1 = ∠ ( ze + zc1 ) = 0.00175158 rad Z c 2 = −3.388518
α 2 = ∠ ( ze + zc 2 ) = arctan
Im( ze ) 2.008100715·10 −3 = arctan Re( ze ) − zc 2 2.473822 + 3.388518
α 2 = ∠ ( ze + zc 2 ) = 0.000467723 rad Zc3 = 51891531.495755
α3 = ∠ ( ze + zc3 ) = arctan
Im( ze ) 2.008100715·10−3 = arctan Re( ze ) − zc3 2.473822 + 51891531.495755
α3 = ∠ ( ze + zc3 ) = 3.8698·10−11 rad
66
POLOS Z p1 = 0.6732
α 4 = ∠ ( z e + z p1 ) = arctan
Im( z e ) 2.008100715·10 −3 = arctan Re( z e ) − z p1 0.904838 − 0.6732
α 4 = ∠ ( z e + z p1 ) = 0.00866892 rad
Z p2 = 1
α 5 = ∠ ( ze + z p 2 ) = arctan
Im( ze ) 2.008100715·10−3 = arctan Re( ze ) − z p 2 0.904838 − 1
α 5 = ∠ ( ze + z p 2 ) = 3.120494 rad El polo que vamos a eliminar, se contrarresta con el cero que añadimos así que podemos prescindir de ellos para el análisis:
α1 + α 2 + α 3 − α 4 − α 5 − α pd = π 0.00175158 + 0.000467723 + 3.8698·10−11 − 0.00866892 − 3.120494 − α pd = π
α pd = −π + 0.00175158 + 0.000467723 + 3.8698·10−11 − 0.00866892 − 3.120494 α pd = −6.268536 α pd
Im( ze ) 2.008100715·10−3 = ∠ ( ze + z pd ) = arctan = arctan Re( ze ) − z pd 0.904838 − z pd
α pd
2.008100715·10−3 = arctan 0.904838 − z pd
( 0.904838 − z ) tan(α pd
pd
) = 2.008100715·10−3
( 0.904838 − z ) = 2.008100715·10 tan(α )
−3
pd
pd
2.008100715·10−3 z pd = − − 0.904838 = 0.767768 tan(−6.268536) Por lo que el controlador PD a añadir quedará del siguiente modo: GPD ( z ) = K d
67
z − 0.9996629 z − 0.767768
El nuevo lugar de las raíces queda de la siguiente manera manera::
Veamos por curiosidad como ha cambiado la respuesta en el tiempo:
Vemos que el sistema se comporta de una forma adecuada pero lenta, además podemos observar un poco de error de aliashing. Vemos que ocurre si añadimos los polos conjugados al lugar comentado comentado.. Observamos que ahora el lugar de las raíces si los contiene, o pasa muy próximo.
En la imagen de la izquierda, tenemos que para situar el polo en 0.905 ± 0.00203 j (muy aproximado al polo que buscamos), necesitamos una ganancia de 113. 113. Si 68
aplicamos esta ganancia obtenemos la respuesta de la derecha. Podemos Podemos ver que no tiene ninguna sobre oscilación oscilación. A continuación podemos ver la respuesta del mismo caso según CC:
Dado que nuestro sistema no es de segundo orden, en los polos encontrados no coinciden del todo los valores de δ y de ωn . Tampoco se llega al valor final de 1 con esta ganancia, por lo que probaremos a aumentar un poco la ganancia buscando un punto de los polos que nos dé un resultado bueno para nuestro el sistema sistema,, con sobre oscilación pequeña pero asegurando que el sistema alcanza el valor buscado buscado:
Probaremos una ganancia entre 148 y 157 por ejemplo 150, la cual deducimos a partir del diagrama anterior que nos dará un resultado magnifico.
69
Con este valor de ganancia Kd = 150 tenemos una sobre oscilación de entre 1 y 2 y hemos superado tanto el amortiguamiento establecido en 0.4 como la frecuencia 25. 25. Al habernos llevado el sistema a frecuencias de 106 rad/s, rad/s, tenemos un sistema muy rápido rápido.. Ve Veamos amos qué respuesta obtenemos:
Podemos ver que hemos satisfecho ambos requisitos ya que no llegamos al 30% de sobre oscilación y el tr < 0.1s . Ahora nos queda comprobar el error en régimen estacionario. Al igual que en el caso continuo, el error discreto equivale a la expresión 1 E ( z ) = R( z ) . Sin embargo en este caso tenemos que el teorema del valor final 1 + G( z ) establece:
e∞ = lim e(k ) = lim (1 − z −1 ) E ( z ) = lim (1 − z −1 ) R( z ) k →∞
70
z →1
z →1
1 1 + G( z )
1 Z {L−1{ f ( s )}} Tz → , que 2 s ( z − 1)2 equivalen a una entrada en rampa, podemos averiguar el error en velocidad del sistema: En consecuencia y conocido el par de transformadas
G6 ( Z ) = K d
z − 0.9996629 −1.062405 ⋅10−13 ( z 3 − 5.189153·107 z 2 − 1.883729· 1.883729·10 108 z − 4.248368·107 ) z − 0.767768 ( z − 1)( z − 0.996629)( z − 0.67032) e∞ = lim(1− z−1 ) z→1
Tz e∞ = lim z→1 (z −1) Tz e∞ = lim z→1 (z −1)
e∞ = lim z→1
Tz 1 Tz 1 = lim (1− z−1 ) 2 2 z → 1 (z −1) 1+G(z) (z−1) 1+G(z)
1 z −0.9996629 −1.062405⋅10 (z −5.189153·10 5.189153·107 z2 −1.883729· 1.883729·10 108 z −4.248368·10 4.248368·107 ) 1+ Kd z −0.767768 (z −1)(z −0.996629)(z −0.67032) −13
1 5.189153·107 z2 −1.883729· 1.883729·10 108 z −4.248368·10 4.248368·107 )) ( z −0.9996629) ( −1.062405⋅10 (z −5.189153·10 −13
1+ Kd
3
3
( z −0.767768) (z −1)(z −0.9996629)(z −0.67032)
( z −0.767768) (z −1) (z −0.9996629)(z −0.67032) Tz (z− z −1) ( z −0.767768) (z −1)(z −0.9996629)(z −0.67032) + Kd ( z −0.9996629) ( −1.062405⋅10−13(z3 −5.189153·10 5.189153·107 z2 −1.883729· 1.883729·10 108 z −4.24836 4.248368·10 4.248368· 8·107 ))
(
Tz( z −0.767768) (z −0.9996629)(z −0.67032)
e∞ = lim z→1
(( z −0.767768) (z −1)(z −0.9996629)(z −0.67032) +K ( z −0.9996629) ( −1.062405⋅10
−13
d
e∞ =
T ( 0.232232) (0.0003371)(0.32968)
( ( z −0.767768) (z −1)(z −0.9996629)(z −0.67032) + K ( 0.0003371) ( −1.062405⋅10
−13
d
e∞ =
e∞ =
T ( 0.232232) (0.0003371)(0.32968)
(
T0.0000258091
(
Kd ( 0.0003371) ( −1.062405⋅10 +5.51298·10 5.51298·10−6 + 2.00127· 2.00127·10 10−5 + 4.51349· 4.51349·10 10−6 ) −13
T 0.0000258091 −13 Kd ( 0.0003371) ( −1.062405⋅10 +5.51298· 5.51298·10 10−6 + 2.00127· 2.00127·10 10−5 + 4.51349· 4.51349·10−6 )
(
e∞ =
T 0.0000258091 T = 2548,744 Kd −1.01262·10−8 Kd
e∞ = 0.001≥
0.001 2.548744 2548,744 = Kd Kd
2.548744 2.548744 →Kd ≥ = 2548.744 Kd 0.001
)
(1−5.189153·10 5.189153·107 −1.883729 1.883729·10 ·10 ·108 − 4.248368·10 4.248368·107 ))
Kd ( 0.0003371) ( −1.062405⋅10−13(1−5.189153·10 5.189153·107 −1.883729·10 1.883729·108 − 4.248368·10 4.248368·107 ))
e∞ =
71
(z3 −55.189153· .189153·10 .189153· 107 z2 −1.883729· 1.883729·10 108 z −4.248368· 4.248368·10 107 ))
)
)
)
)
)
Como podemos ver en la siguiente imagen, no podemos obtener el error que queremos dentro del margen, sin sacrificar la sobre oscilación. oscilación. Por tanto no podemos obtenerlo mediante el incremento del K d . También cómo podemos ver, el erro errorr depende del periodo de muestreo por lo que podríamos haber obtenido un menor error estacionario, con un periodo de muestreo menor.
La solución que usaremos será, será usar una ganancia de lo más alta que podamos y por otro lado implementa implementar un controlador PI, que nos permitirá mejorar el régimen permanente sin alterar el estacionario. De la imagen anterior, observamos que podemos aplicar una ganancia de 320 sin obtener un sistema que sobre oscile demasiado. Veamos el error que nos quedaría por corregir: e∞ =
2.548744 = 0.0079648 320
Por tanto nuestro error en estado estacionario de velocidad, aún supera unas 8 veces el error máximo que queremos tener. El sistema en con esta ganancia nos da la siguiente respuesta en el tiempo:
72
Podemos ver que seguimos teniendo una respuesta muy buena, sobre todo cabe destacar la rapidez que se ha obtenido, en 0.02 segundos ya se ha alcanzado el punto de máxima sobre oscilación. oscilación Vamos ahora con el diseño del controlador integral, añadiremos un controlado que cumpla cumpla:: − zci < − z pi z − zci Gci ( z ) = , siendo z − z pi zci > z pi Debemos escoger los valores para el polo y para el cero, muy próximos a 1, 1, de forma, que solo modificarán las bajas frecuencias (al igual que en el diagrama de Bode) manteniendo invariable el régimen transitorio del sistema. e∞ = e∞ =
T 0.001 → 0.001 ≥ → K ve ≥ 1 K ve K ve
T 2.548744 2548, 744 = = 0.00796483 Kd 320 Kv =
0.001 = 0.125552 0.00796483
necesitamos un K v ≥ 1 y tenemos K v = 0.125552 K ve 1 = ≈8 K v 0.125552
Usaremos la relación entre el valor de K p que necesitamos tener y el que tenemos actualmente para calcular el cero del controlador, a partir de un polo que tomaremos cercano a 1 por ejemplo: z pi = −0.99 en los casos de una relación mayor entre las constantes de error, necesitaremos comen comenzar, ar, más cerca de la circunferencia unidad.
1 + zci 1 + zci =8→ = 8 → zci = 8·0.01 + 1 → zci = 0.08 + 1 = 1.08 1 + z pi 1 + 0.99 1 + 1.08 I= 1 + 0.99
73
El controlado controlado Integral resultante parece bueno en principio, probaremos a ver qué resultado obtenemos:
Podemos ver que el resultado es bastante satisfactorio. Cabe comentar que podríamos haber colocado el polo y cero en el lado positivo del eje, sin embarg embargo o hubiesen interferido con los polos importantes del sistema y desajustado lo que habíamos conseguido con el PD PD.. Habiendo colocado estos dos en la parte negativa del eje real del sistema, se cancelan uno a otro, y además no interfieren en el comportamiento del sistema, ya que no modifican el lugar de las raíces de manera significativa significativa.. Como los polos se encuentran en el lugar que queremos para una ganancia de 1, el polo nuevo no se moverá prácticamente, pero intervendría a partir de ahora si se quisiese aum aumentar entar la ganancia del sistema, ya que por encima de 227 el sistema se volvería inestable. Esto podemos compro arlo en las imágenes siguientes:
74
Podemos ver que seguimos cumpliendo con la sobre oscilación y con el tiempo de subida, por lo que solo nos queda asegurarnos de que cumplimos el error en estado estacionario. G6 ( Z ) =
z + 1.08 z − 0.9996629 −1.062405 ⋅10−13 ( z 3 − 5.189153·107 z 2 − 1.883729·10 1.883729·108 z − 4.248368·107 ) Kd z + 0.99 z − 0.767768 ( z − 1)( z − 0.996629)( z − 0.67032)
e∞ = lim(1− z−1 ) z→1
Tz e∞ = lim z→1 (z −1) Tz e∞ = lim z→1 (z −1)
e∞ =lim z→1
Tz 1 Tz 1 = lim (1− z−1 ) 2 (z −1)2 1+G(z) z→1 1 + G (z) (z−1)
1 z +1.08 z −0.9996629 −1.062405⋅10 (z −5.189153· 5.189153·10 107 z2 −1.883 1.883729· 729·10 729· 108 z −4.248368·10 4.248368·107 ) 1+ Kd z +0.99 z −0.767768 (z −1)(z −0.996629)(z −0.67032) −13
3
1 5.189153·10 107 z2 −1.883729· 1.883729·10 108 z −4.248368·10 4.248368·107 )) ( z +1.08)( z −0.9996629) ( −1.062405⋅10 (z3 −5.189153· −13
1+ Kd
( z +0.99)( z −0.767768) (z −1)(z −0.9996629)(z −0.67032)
( z+0.99)( z−0.767768) (z−1)1)(z−0.9996629)(z−0.67032) Tz (z− z −1) ( z −0.767768) (z −1)( 1)(z −0.9996629)(z −0.67032)+Kd ( z +1.08)( 0.9996629 6629) ( −1.062405⋅10−13(z3 −5.189153·107 z2 −1.883729·108 z −4.2 4.248368· 48368·107)) ) z−0.999
(
Tz( z +0.99)(( z −0.767768) (z −0.9996629)(z −0.67032) e∞ =lim z→1 Kd ( z −0.767768) (z −1)( 1)(z −0.9996629)(z −0.67032)+Kd ( z +1.08)( z −0.9996 0.9996629 629) ( −1.062405⋅10−13(z3 −5.189153·107 z2 −1.883729·108 z −4.248368· 4.248368·107))
(
e∞ =
T ( z+ z +0.99)( ) 0.232232) (0.0003371)(0.32968)
( ( z −0.767768) (z −1)(z −0.9996629)(z −0.67032) +K ( z +1.08)() 0.0003371) ( −1.062405⋅10
−13
d
e∞ =
Kd (2.08)( 0.0003371) ( −1.062405⋅10−13(1−5.189153·107 −1.883729· 1.883729·108 −4.248368·107 ))
e∞ =
e∞ =
T0.000000258091
(
Kd (2.08)( 0.0003371) ( −1.062405⋅10−13 +5.51298·10−6 +2.001 2.00127· 27·10−5 +4.51349·10−6 )
T 0.000000258091 Kd (2.08)( 0.0003371) ( −1.062405⋅10−13 +5.51298· 5.51298·10 10−6 + 2.00127· 2.00127·10−5 + 4.51349·10−6 ) 2.00127·10
(
e∞ = e∞ =
75
)
(1−5.189153·107 −1.88 1.883729· 3729·108 −4.248368·107 ))
T ( 0.01)( 0.232232) (0.0003371)(0.32968)
(
)
T 0.000000258091 T = 12.253579 Kd 2.10625· 2.10625·10 10−8 Kd
0.001 0.012253579 0.012253579 12.253579 = = = 0.0000382924 Kd Kd 320
)
)
)
)
Podemos ver que cumplimos perfectamente con el error máximo que teníamos para el sistema. Por lo que hemos finalizado el diseño diseño. Nuestro uestro sistema final es: G9 ( Z ) =
z + 1.08 320 ( z − 0.9996629 ) −1.062405 ⋅10−13 ( z 3 − 5.189153·107 z 2 − 1.883729·108 z − 4.2 4.248368·10 48368·107 ) · z + 0.99 z − 0.767768 ( z − 1)( z − 0.996629)( z − 0.67032) Gc ( Z ) =
z + 1.08 320 ( z − 0.9996629 ) · z + 0.99 z − 0.767768
Por tanto, tenemos el siguiente controlador discreto para nuestro sistema, que nos pe permite rmite cumplir con los requerimientos satisfactoriamente satisfactoriamente:
θ R (Z ) E(Z )
+
−θ (Z ) S
U (Z )
GcPID ( z )
0.001081191K a 1 − eTs U (s) s ( 2.966004s + 1)(( 0.0025s + 1) s Bloqueador de orden cero A/D D
K a = 2500000 320 ( z − 0.9996629 ) z + 1.08 320 ( z − 0.9996629 ) GcPID ( Z ) = · PD = z + 0.99 z − 0.767768 z − 0.767768 z + 1.08 I= z + 0.99
76
θ S (s)
9. Análisis de los sistemas obtenidos obtenidos. En este apartado, realizaremos un análisis de los sistemas que hemos obtenido y que cumplen con los requerimientos impuestos. Estudiaremos la señal de error y de control de cada uno de ellos, ellos y veremos las implicaciones de construir dos de ellos ellos.. Concretamente analizaremos, el controlador PD diseñado mediante la aproximación a Concretamente sistema de segundo orden, el PID diseñado a través de su respuesta en frecuencia, y por último el PID discreto, diseñado mediante el lugar de las raíces en el plano Z.
R( s)
+
E ( s)
Gr ( s )
U (s)
Y (s)
Gb ( s )
G (s)
−
Como se ha demostrado en ocasiones anteriores, sabemos que la señal de error es R( s ) E (s) = donde R(s) es la señal de referencia y G(s) es la función de transferencia 1 + G( s) completa del sistema. También sabemos qué U (s) = E (s)Gr (s) y por tanto tenemos:
Ger ( s) =
R( s)Gr ( s) R( s) U ( s) = 1 + G( s) 1 + G( s)
Estas dos señales, son las señales de error y de control. Ellas nos pueden dar información interesante sobre el sistema. De ellas podemos extraer el esfuerzo que el sistema ha de re realizar alizar para cumplir los requisitos que le hemos impuesto, y en base a ese esfuerzo podemos saber si el voltaje necesario para controlar el sistema es admisible o no por la electrónica del circuito. Para realizar estos análisis, se empleará una señal de ref referencia, erencia, escalón unitario, y se analizará la función resultante de las ecuaciones anteriores, mediante la herramienta CC, utilizando un impulso en lazo abierto. abierto Dee esta manera pode podemos mos ver la evolución en el tiempo, de las funciones de transferencia equival equivalentes entes a las señales que queremos estudiar estudiar,, ante un impulso unitario. unitario.
77
9.1 PD En primer lugar comenzaremos analizando, el controlador PD obtenido. Recordemos que en nuestro caso el amplificador K a estaba ya integrado en el circuito original, original, no obstante, forma parte del regulador del sistema pero no lo tendremos en cuenta a la hora de ver la señal de control, supondremos sencillamente que esta señal será amplificada posteriormente posteriormente:
K a = 2500000 GcPD ( s ) = 2500000 (1 + 0.046292 ) G pd ( s ) = 1 + 0.046292s El sistema regulado queda como sigue:
Vamos las señales de control y error del sistema:
Ger ( s) =
R( s)Gr ( s) R( s) U ( s) = 1 + G( s) 1 + G( s)
A continuación, vemos la señal de error en rojo y la señal de control en azul. Podemos observar, que los valores que alcanzan ambas señales no son valores, desorbitados, no llegan a los 5 voltios en ningún momento, por lo que en principio son valores aceptables para la construcción del sistema.
78
9.2 PID (Diseñado por respuesta en frecuencia) Ahora, analizaremos el controlador PID diseñado mediante el sistema de la respuesta en frecuencia y diagramas de Bode. Este controlador se corresponde con una red de adelanto y otra red de atraso de frecuencias, y una amplificación llevada a cabo por el am amplificador plificador operacional, empleado en la planta K a . Por tanto partimos de las siguientes funciones del sistema: K a = 2000000 1 + 3.0227553s 1 + 0.21144s 1 + 0.21144s GcPID ( s ) = 4.8 0.099415 Gav ( s ) = 0.099415 1 + 14.509225s 1 + 0.0210199s 1 + 0.0210199s 1 + 3.0227553s Gret ( s ) = 4.8 1 + 14.509225s
El sistema completo queda como sigue: G5 ( s ) = 4.8
1 + 3.0227553s 1 + 0.21144s 2162.382 0.099415 1 + 14.509225s 1 + 0.0210199s s ( 2.966004 s + 1) ( 0.0025s + 1)
De aquí extraemos que tenemos tenemos que diseñar una red de adelanto y otra de atraso atraso.. Estas redes nos servirán como reguladores para el sistema. El esquema general de una red de adelanto o atraso es el siguiente: C2
R2 R4 4
C1
R1
R3
E1 ( s )
VCC
E2 (s)
E3 ( s )
Dado este sistema, y dependiendo del valor dad dado o a α tendremos una red de adelanto o una red de atraso. A continuación se pasa a explicar las distintas ecuaciones que definen estos valores, analizaremos la red anterior como un controlador: Gc ( s ) =
79
E3 ( s ) E3 ( s ) E2 ( s ) R2 · R4 1 + R1C1s 1 + Ts = = , y sabiendo que: Gc ( s ) = K cα E1 ( s ) E2 ( s ) E1 ( s ) R1· R3 1 + R2C2 s 1 + α Ts
Esto nos lleva a las siguientes ecuaciones: α = α T = R2 ·C2 T = R1C1 T R1·C1 Si R1·C1 > R2 ·C2 ⇒ α < 1 → R ed de adelanto α T = T2C2 → Si R1·C1 < R2 ·C2 ⇒ α > 1 → R ed de atraso K α R C c = 4 1 R2 · R4 K c = K cα = α R3C2 R1·R3 En el caso de nuestra red, la ganancia es 1 en los dos casos ya que no hemos tenido que añadir ninguna ganancia extra. Por tanto comenzamos con la búsqueda de los valores, comenzamos con la red de adelanto:
K c = 1 → R4 ·C1 = R3·C2
α = 0.099415 =
R2 ·C2 R2 · R4 = R1·C1 R1·R3
T = R1·C1
α T = R2·C2 Vamos a suponer R2 = 10K Ω y por tanto tenemos: 0.0210199 α T = 10000·C2 = 0.0210199 → C2 = = 2.10199·10−6 = 2.1µ F 10000 Podríamos usar un capacitor de voltaje máximo de 7960 con capacidad de 2.093 µ F que está bastante cerca del valor necesario. De modo que αT = 2.093·10−6·10000 K c = 1 → R4 ·C1 = R3 ·0.0210199 → R3 =
1 = 477783.0865Ω ≈ 470 K Ω = 470000Ω 2.093·10−6
Podemos ver que la resistencia R3 puede construirse, por encima o por debajo de su valor en una pequeña franja. Debido a que de esta resistencia depende la ganancia, la ponemos un poco más baja y de esta manera ganaremos un poco de ganancia en el peor caso. Para continuar, estableceremos por ejemplo la resistencia R1 = 1M Ω , continuamos por tanto: R3 = 470000Ω R2 = 10 K Ω R1 = 1M Ω C2 = 2.093µ F
α T = 0.02093 T = R1 ·C1 = 1M Ω·C1 = 0.21144 → C1 =
0.21144 = 0.2114µ F ≈ 0.271µΩ 1M Ω
10000·R4 10000·2.093·10−6 = = 0.0772325 1M Ω·0.271µ F 1M Ω·470000 772325·470000 R4 = = 36299275 ≈ 4·9.1M Ω 10000 R4 = 4·9.1M Ω = 36400000Ω
α = 0.099415 =
80
Recopilemos los resultados: R3 = 470000Ω R2 = 10 K Ω R1 = 1M Ω C2 = 2.093µ F R4 = 36.4 M Ω C1 = 0.271µΩ De modo que nuestro controlador queda: Gc ( s ) = Gc ( s ) =
R2 ·R4 1 + R1C1s 1 + Ts , que equivale a: Gc ( s ) = K cα R1 ·R3 1 + R2C2 s 1 + α Ts
10 K Ω·36.4M Ω 1 + (1M Ω·0.271µ F ) s 1 + 0.271s · → 0.0774468 1M Ω·4.7 M Ω 1 + (10 K Ω·2.093µ F ) s 1 + 0.02093s
Como podemos ver, el sistema real no se corresponde del todo con los valores del controlador teórico debido a los valores de mercado a los que tenemos que ajustarnos. Como podemos ver el valor de α nos queda un poco bajo, podemos intentar corregirlo un poco, modificaremos el valor de R1 = 100K Ω y así nos permitirá reajustar algunos valores valores:: R3 = 470000Ω R2 = 10 K Ω R1 = 1M Ω C2 = 2.093µ F
α T = 0.02093 T = R1 ·C1 = 1M Ω·C1 = 0.21144 → C1 =
0.21144 = 2.114µ F ≈ 2.167 µΩ 100 K Ω
10000·R4 10000Ω·2.093·10−6 F = = 0.0965851 −6 100 K Ω·2.167·10 F 100 K Ω·470000 9658.51·470000 R4 = = 453949.97 ≈ 390 K Ω + 68 K Ω = 458 K Ω 10000
α = 0.099415 =
Por lo que con estos cálculos, tenemos: R3 = 470000Ω R2 = 10 K Ω R1 = 1M Ω C2 = 2.093µ F R4 = 470 K Ω C1 = 2.167 µΩ
Gc ( s ) = Gc ( s ) =
R2 ·R4 1 + R1C1s 1 + Ts , que equivale a: Gc ( s ) = K cα R1 ·R3 1 + R2C2 s 1 + α Ts
10 K Ω·458 K Ω 1 + (100 K Ω·2.167 µ F ) s 1 + 0.2167 s · → 0.0974468 100 K Ω·470 K Ω 1 + (10 K Ω·2.093µ F ) s 1 + 0.02093s
Gc ideal ( s ) = 0.099415
1 + 0.21144s 1 + 0.0210199s
Gc real ( s ) = 0.0974468
1 + 0.2167 s 1 + 0.02093s
Podemos ver, que no es exactamente la misma red la que hemos conseguido, pero sí que se acerca mucho. Más adelante veremos el resultado de estas redes. 81
Vamos a continuación con la red de atraso de fase, comenzamos igualmente por comprobar los valores que queremos conseguir, y estableciendo R2 = 10K Ω
K c = 1 → R4 ·C1 = R3 ·C2
α = 4.8 =
R2 ·C2 R2 ·R4 = R1 ·C1 R1 ·R3
T = R1 ·C1 = 3.0227553
α T = R2 ·C2 = 14.509225 Hallemos los valores de los condensadores y resistencias en función del valor de R2 :
K c = 1 → R4 ·C1 = R3 ·C2
α = 4.8 =
10 K Ω·C2 10 K Ω·R4 = R1 ·C1 R1 ·R3
T = R1 ·C1 = 3.0227553
α T = 10 K Ω·C2 = 14.509225 → C2 =
14.509225 = 1450.809µ F ≈ 14·103.715µ F = 1452.01µ F 10000
Necesitaríamos 14 condensadores en paralelo de tipo 1452.01 µ F construcción del condensador C2 . Hallemos ahora el valor de R3 : K c = 1 → R4 ·C1 = R3 ·1452.01µ F → R3 =
1 = 688.700491Ω ≈ 680Ω 1452.01µ F
Para hallar el resto de valores, establecemos el valor de R1 = 1M Ω : 3.0227553 T = 1M Ω·C1 = 3.0227553 → C1 = = 3.0227 µ F ≈ 3.066µ F 1M Ω Por tanto, veamos el valor de α y hallemos en funci función ón de él, él, el valor de R4 : α=
10 K Ω·R4 10 K Ω·1452.01µ F = 4.7358 → = 4.7358 1M Ω·3.066 µ F 1M Ω·680Ω 1M Ω·680Ω 4.7358 R4 = = 3220340.4Ω ≈ 330 K Ω 10 K Ω
En resumen tenemos: R3 = 680Ω R2 = 10 K Ω R1 = 1M Ω C2 = 1452.01µ F R4 = 330 K Ω C1 = 3.066 µΩ Por lo que finalmente la red de retraso final quedaría del siguiente modo:
82
para la
Gc ( s ) = Gc ( s ) =
R2 ·R4 1 + R1C1s 1 + Ts , que equivale a: Gc ( s ) = K cα R1 ·R3 1 + R2C2 s 1 + α Ts
10 K Ω·330 K Ω 1 + (1M Ω·3.066 µ F ) s 1 + 3.066s · → 4.7358 1M Ω·33K Ω 1 + (10 K Ω·1452.01µ F ) s 1 + 14.5201s
Por tanto, el controlador queda: Gc real ( s ) = 4.7358
1 + 3.066 s 1 + 0.2167 s 0.0974468 1 + 14.5201s 1 + 0.02093s
Veamos su respuesta en el tiempo:
Podemos ver que tenemos un comportamiento del sistema, más que aceptable, faltaría comprobar el error en régimen estacionario, pero los requisitos de sobre oscilación y tiempo de subida se cumplen sin problemas. Veamos la diferencia en señal de error y de control entre el controlador real en rojo y el ideal en azul:
Podemos observar que las señales señales de error y control obtenidas son muy parecidas. De hecho mantienen valores satisfactorios, porque no llegan a alcanzar esfuerzos desproporcionados, por tanto podemos concluir que el sistema parece factible para su construcción.
83
9.1 P PID ID Discreto (Diseñ (Diseñado ado por lugar de las raíces raíces) Como se comentó anteriormente, un sistema discreto, se realiza siguiendo un esquema como el siguiente:
θ R (Z ) E(Z )
+
−θ (Z ) S
U (Z )
GcPID ( z )
0.001081191K a 1 − eTs U (s) s ( 2.966004s + 1)(( 0.0025s + 1) s
θ S (s)
Bloqueador de orden cero A/D D
Para su construcción, requeriremos un PC, o bien una placa con un micro controlador. En el ca caso so de este tipo de sistemas, los factores de mayor problemática, están en la frecuencia de reloj, y en la capacidad del microprocesador para ejecutar las instrucciones suficientemente rápido para obtener las señales de control en función de su necesidad. K a = 2500000 320 ( z − 0.9996629 ) z + 1.08 320 ( z − 0.9996629 ) GcPID ( Z ) = ⋅ PD = z + 0.99 z − 0.767768 z − 0.767768 z + 1.08 I= z + 0.99
Para poder implementar el PID como un controlador digital, el proceso a seguir, sería; en primer lugar transformar la función de control al dominio del tiempo, en segundo lugar implementar un algoritmo que ejecute la parte integral y derivada de esa función y las emplee para corregir el error de la salida del sistema respecto a la entrada del mismo, y a continuación devuelva a la planta la señal de control construida. Por último cabe comentar, que habrá que ejecutar este algoritmo, algoritmo, en un sistema, que utilice una función lanzada de forma periódica, cada ( Ts segundos) y que ejecute el código antes mostrado. Es importante comentar que no debe exceder su ejecución, el tiempo de muestreo, pues si esto ocurri ocurriese, ese, se solaparían ejecuciones diferentes del bucle y el controlador podría provocar resultados imprevistos. Por falta de tiempo, este método no se aplicará, pero si comentaremos algunos aspectos importantes de la construcción, derivados de los comentario comentarioss anteriores: Se ha usado un periodo de muestreo para el sistema Ts = 0.001s , esto implicará que el reloj mínimo debería ser capaz de lanzar una función periódica que ejecute la lógica de nuestro regulador con una frecuencia de 1Mhz. Si suponemos que un algoritmo de control, tiene una media de unas 50 instrucciones en ensamblador y que usamos un micro que es capaz de ejecutar una instrucción ccompleta ompleta por ciclo (como suele ocurrir en procesadores segmentados cuando el tiempo tiende a infinito). Necesitaríamos un 84
controlador con una velocidad de 50Mhz como mínimo para poder ejecutar nuestro algoritmo correctamente. Por lo que la implementación di digital gital del controlador no sería muy sencilla ya que la mayoría de los micro microcontroladores controladores suelen trabajar con relojes del orden de 1010-20Mhz. 20Mhz. Necesitaríamos un microcontrolador potente un PC o una FPGA para poder implementar el controlador digital necesario. Según se vio al realizar el diseño digital del sistema, podríamos haber usado una periodo de muestreo mayor, pero con 0.005 segundos, segundos la señal original del sistema ya no se podría construir. Por tanto el mecanismo óptimo, óptimo, sería utilizar el teorema de Nyqui Nyquist st para poder determinar la frecuencia mínima de muestreo y así optimizar el reloj necesario para ejecutar el regulador regulador. Veremos a continuación la señal de error en rojo y de control en azul producidas por el sistema:
Como podemos ver en la imagen, lo loss valores de control son enormes frente a los valores de error. Esto seguramente será debido a la discretización de los valores, que hace que las decisiones de control a tomar hayan de ser mucho más drásticas. Veamos las señales por separado:
Como conclusión podemos extraer que la construcción de nuestro sistema de forma digital, conlleva más complicación que la construcción analógica analógica.. No obstante cabe 85
destacar las ventajas que también tenemos con el control digital. digital. Para empezar tenemos un menor menor coste en trabajo, trabajo, ya que los sistemas están más probados (si se usan placas entrenadoras con controladores o bien PCs) son sistemas con una baja pro probabilidad babilidad de error error.. Por tanto los errores se reducen a fallos de programación. El El coste sería la programa programación ción y por tanto como materia prima solo necesitamos horas de trabajo, además los algoritmos también están estandarizados y sólo sería necesario ajustar los valores de nuestro regulador. Un n error de código no implicaría rehacer todo el trabajo, mientras que que el hecho de tener que construir un sistema físico conlleva los inconvenientes de que una vez fabricado no es fácil la vuelta atrás. Un sistema analógico puede puede que no reaccione del todo como estaba planeado por problemas de tolerancias de los componentes, u otros aspectos que nos llevarían a tener que construir de nuevo el controlador controlador..
86
10. ANEXO A Valores comerciales de resistencias
Calculo de resistencias
87
11. ANEXO B Valores comerciales de condensadores
88
