CONTROL PID DE UN SERVOMECANISMO NO LINEAL USANDO LA FUNCION DESCRIPTIVA Y LA BIFURCACION DE BOGDANOV-TAKENS Manuel F. Pérez Polo, Manuel Pérez Molina, Javier Gil Chica, José Ángel Berná Galiano Departamento de Física Ingeniería de Sistemas y Teoría de la Señal. Universidad de Alicante. Escuela Politécnica Superior. Apartado 99, E-03080. Alicante. E-mail: [email protected]; ; [email protected] ; [email protected] ; [email protected]

complejos conjugados puros) es una excepción, ya que ha sido usada en el análisis y diseño de muchos sistemas de control [11], [6] sin embargo, la bifurcación de Bogdanov-Takens (cuando la matriz de la parte lineal del sistema tiene dos autovalores cero) se ha usado mucho menos que la anterior [11], [6], [10], [1], [2], [12]. En este trabajo se presenta un método de control de un servomecanismo no lineal de posición [8] basado en un regulador PID cuyos parámetros se ajustan teniendo en cuenta la bifurcación de Bogdanov-Takens [3], [5], [7], [13]. Para ello se obtiene el modelo del sistema reducido a la variedad central, a partir del cual se estudia la estabilidad usando el segundo método de Lyapunov, lo que permite determinar la constante de acción proporcional e integral del regulador PID. La acción derivada del regulador PID permanece desconocida, por lo que se ha utilizado la función descriptiva, lo que permite estimar tanto la frecuencia de auto-oscilación del sistema como la constante de acción derivada del regulador PID [8], [13] Una vez alcanzado el régimen de autooscilación, se comprueba que una variación armónica de la entrada de referencia puede dar lugar a comportamiento caótico [9], [10], [11], que se ha usado para mejorar el funcionamiento del regulador PID. Finalmente, la robustez del diseño se ha comprobado suponiendo que existe ruido aleatorio en el par de perturbación aplicado en el eje de salida.

Resumen En este trabajo se investiga como ajustar los parámetros de un regulador PID, que controla la posición del eje de un servomecanismo con una no linealidad cúbica. Se demuestra a partir de las ecuaciones del modelo y del regulador PID, que en el jacobiano del sistema linealizado pueden aparecer dos ceros, lo que da lugar a una bifurcación de Bogdanov-Takens, cuyas condiciones de estabilidad se estudian con el teorema de la variedad central y el segundo método de Lyapunov, lo que permite determinar la acción proporcional e integral del regulador. La acción derivada se calcula a través de la función descriptiva de la no linealidad cúbica, lo que permite obtener una aproximación de la frecuencia de auto-oscilación. Se supone que la entrada de referencia se puede variar senoidalmente y que existe una pequeña perturbación armónica en el eje de salida del servomecanismo. En este contexto, pueden aparecer oscilaciones caóticas que se utilizan para obtener señales de control pequeñas. El ajuste del regulador PID se comprueba en presencia de ruido aleatorio. Palabras clave: Servomecanismo, no linealidad cúbica, ajuste regulador PID, bifurcación de Bogdanov-Takens, función descriptiva, caos

1 INTRODUCCIÓN

2 MODELO DEL SISTEMA

El control de dispositivos electromecánicos se ha centrado casi exclusivamente en la utilización de modelos lineales, a pesar de que es conocido que en tales dispositivos aparecen fuertes no linealidades cuando funcionan en régimen de grandes desviaciones respecto al estado de equilibrio [4]. Una posible solución a este problema es usar la teoría de bifurcaciones, a pesar de que en sistemas de dimensiones mayores que tres, se ha usado poco en el análisis y diseño de un posible sistema de control. La bifurcación de Hopf (que aparece cuando en la matriz de la parte lineal del sistema tiene un par de autovalores

El sistema está formado por un cuerpo, cuya posición se desea controlar, conectado a un tren de engranajes considerado ideal con relación de transmisión n = N1/N2. El tren de engranajes está conectado a un motor DC controlado por armadura, al cual se le aplica el voltaje de entrada ua(t) a través de un amplificador de constante Ka. El sistema se realimentado con un potenciómetro de constante Ks cuyo voltaje de salida se aplica al regulador PID conectado a una no linealidad cúbica de constante Kne. Los

1

u2 ( t ) = K ne u1 ( t )

valores de los parámetros se indican en la leyenda de la figura 1. A partir de la figura 1.las ecuaciones del modelo matemático del servomecanismo y sistema de control son las siguientes: Potenciómetro: u ( t ) = K s [θ r ( t ) − θ ( t )]

3

siendo Jm y JL los momentos de inercia del motor y de la carga respectivamente. Bm y BL son la fricción viscosa del motor y de la carga y Jme y Bme son el momento de inercia y la fricción viscosa reducidos al eje del motor. Los valores, Tm(t), T’m(t), T’L(t) and Td(t) representan el par motor, el par motor para mover el eje más la fricción viscosa, el par de carga y el par de perturbación respectivamente. Par motor asumiendo régimen lineal:

(1)

Regulador PID:

⎡

1

⎣

τi

u1 ( t ) = K p ⎢ u ( t ) +

t

∫

u ( σ ) dσ + τ d

0

du ( t ) ⎤ dt

⎥ ⎦

(2) siendo Kp, τi y τd las constantes del regulador.

Tm ( t ) = K m ia ( t )

Regulador PID (Kp, τi, τd)

Potenciometro K2 + +

+

u1(t)

u(t)

θr(t) = θr1 + θr2sin(ωrt)

θ(t)

(6)

en donde Km es la constante del par.

Elemento no lineal:

+

(3)

Elemento no lineal Kne

n(t) = fna[ X – 0.5] u2(t) θ (t) u (t) m a Tren de engranajes Motor DC Amplificador más carga Ra, La, ub(t) Ka n = N1/N2 Jm, Bm JL, BL

Td(t) = Td2sin(ωdt) Figura 1: Diagrama de bloques del servomecanismo con el sistema de control. Jm = 1.36.10-2 kg.m2, Bm = 9.64.10-3 N.m.s, Ra = 5 Ω, La = 0.01 H, par motor Km; constante de fuerza contra-electromotriz Kb = Km Parametros de la carga: JL = 1.36 kg.m2, BL = 0.136 N.m.s; n = 0.1; Ks = 10/π V/rad y Ka = 10-100; Kne 0.1 V-2. Intensidad de corriente en la armadura: di ( t ) La a + Ra ia ( t ) = ua ( t ) − ub ( t ) dt ub ( t ) = K bω m ( t ) = K b

dθ m ( t )

(5) La entrada de referencia θr(t) y el par de perturbación se definen a través de las ecuaciones:

(4)

θ r ( t ) = θ r 1 + θ r 2 sin ( ω r t ) ; Td ( t ) = Td 2 sin ( ω d t )

dt siendo ub(t) la fuerza contra-electromotriz, θm(t) es el ángulo girado en el eje del motor y is the ωm(t) es la velocidad angular del eje del motor.

(7) en donde los parámetros de las ecuaciones (7) pueden seleccionarse arbitrariamente. Puesto que θm = (N1/N2)θ, y teniendo en cuenta que es posible eliminar la corriente de la armadura ia(t) entre las ecuaciones (4) y (5).y si se tomen derivadas de hasta tercer orden en la ecuación (1), teniendo en cuenta las ecuaciones (5), (6) y

Ecuaciones del par motor: J me

d 2θ m ( t ) dt 2

+ Bme

dθ m ( t ) dt

= Tm ( t ) ; J me = J m + ( N1 N 2 ) J L 2

Bme = Bm + ( N1 N 2 ) BL ; Tm ( t ) = Tm′ ( t ) + TL′ ( t ) + Td ( t ) 2

2

(7) se puede obtener el modelo del sistema de la siguiente forma:

Se puede comprobar que los autovalores de la matriz A (ecuación (11)) son:

d 3u (t ) 2ξ d 2u (t ) 1 du (t ) +τ + = dt 3 dt 2 τ 2 dt ⎡ ⎤ − K2 u13 (t ) + nK s ⎢⎢ Ra Td 2 sin ⎛⎜ωd t ⎞⎟ + Td 2ωd cos ⎛⎜ωd t ⎞⎟ ⎥⎥ J L ⎝ ⎠ ⎝ ⎠ τ me ⎣ a ⎦ 2 ξ K K 2 s s θ ω sin ω t − + θ ω cos ω t −

τ2

( r)

τ r2 r − K sθr 2ωr2 cos (ωrt )

r2 r

λI − A = λ

( r)

J me La ; K ′ = nKm Ka τ 2 Km Kb + Ra Bme J me La

⎛ ⎞ J me La 2ξτ = ⎜⎜ Bme + Ra ⎟⎟ ; K = Kne Ks K ′ L ⎜ J me a ⎟⎠ Km Kb + Ra Bme ⎝ ⎞

ξ = 1 ⎜⎜ Bme + Ra ⎟⎟

J me La 2 ⎜⎝ J me La ⎟⎠ Km Kb + Ra Bme ⎤

⎣

⎦

⎛ ⎞ ⎛ ⎞ ⎢⎜ B ⎟ ⎜ R ⎟⎥ ⎢⎜ me ⎟ + ⎜ a ⎟ ⎥ − R B Km = La J me 4ξ 2 ⎢⎢⎜⎝ J me ⎟⎠ ⎜⎝ La ⎟⎠⎥⎥ a me

(9) donde Km se calcula de forma que 0 < ξ < 1. Se puede comprobar que ξ es adimensional mientras que las unidades de τ, K’ y K son (s), V-1s-1 y V-2s-1 respectivamente. Las ecuaciones (8) y (9) se usarán para analizar la presencia de la bifurcación de Bogdanov-Takens [3], [8], [9], [13].

(t ) =

du ( t )

d u (t )

3

dt

dt

2

; x4

(t ) = u (t ) 1

las ecuaciones (8) se pueden escribir de la forma: 1

0

0

1

−1 τ Kp

2

−2ξ τ K pτ d

2

)] = 0 ..(12)

τ ; for 0 < ξ < 1

2

−α K

τ

p

α +ω 2

2

i

+K

p

i

K

0

p

1

0

ω

−α

0

−2αβ

α −ω

0

b

a

(1 − τ α ) d

;b =

−ω K

2

p

τ

α +ω 2

2

i

2

⎤ ⎡x⎤ ⎥ ⎢ y⎥ ⎥⎢ ⎥ ⎥⎢z⎥ ⎥⎦ ⎢⎣ w ⎥⎦

+Kτ ω p

d

(13) Paso 2. A partir de la s ecuaciones (13) las ecuaciones (11) se pueden escribir de la forma:

(10)

⎡ x1 ( t ) ⎤ ⎡ 0 ⎢ x ( t ) ⎥ ⎢ 0 ⎢ 2 ⎥=⎢ ⎢ x3 ( t ) ⎥ ⎢ 0 ⎢ x ( t ) ⎥ ⎢ K τ ⎣ 4 ⎦ ⎣ p i

2

4

2

; x (t ) =

τ

3

Suponiendo que no hay par de perturbación ( Td2 = 0) y que la entrada de referencia es un escalón (θr2 = 0), introduciendo las variables de estado: ; x2

⎡ x ⎤ ⎡0 ⎢ x ⎥ ⎢0 ⇒⎢ ⎥=⎢ ⎢ x ⎥ ⎢0 ⎢ x ⎥ ⎢1 ⎣ ⎦ ⎣

x = Px

3 ANALISIS DE LA BIFURCACION BOGDANOV-TAKENS

(t ) = u (t )

1−ξ

1

a=

x1

(

+ ( 2ξ τ ) λ + 1 τ

y por consiguiente aparece un doble cero en el punto de equilibrio de las ecuaciones (11), que es el origen. El análisis de la bifurcación se lleva a cabo a través de una serie de pasos para obtener un sistema simplificado reducido a la variedad central [3], [5], [13] el cual mantiene todas las propiedades dinámicas del sistema en un entorno del punto de equilibrio. Paso 1. Sea P la matriz de paso que transforma la parte lineal de la ecuación (11) en su forma canónica de Jordan. Las columnas de P son los vectore propios de los autovalores de la ecuación (12):

(8)

⎡

2

α =ξ τ ; ω =

siendo:

⎛

[λ

λ1, 2 = 0 ; λ3 , 4 = −α ± ω j

du (t ) d 2u (t ) du1 (t ) K p = τ u (t ) + K p + K pτ d dt dt dt 2 i

τ2=

2

0 ⎤ ⎡ x1 ( t ) ⎤ 0 ⎥ ⎢ x2 ( t ) ⎥

⎥⎢ ⎥+ 0 ⎥ ⎢ x3 ( t ) ⎥ ⎥⎢ ⎥ 0 ⎦ ⎣ x4 ( t ) ⎦

⎡ 0 ⎤ ⎢ 0 ⎥ ⎥ x43 ( t ) ⇒ x ( t ) = Ax ( t ) + Bx43 +⎢ 2 ⎢−K τ ⎥ ⎢ 0 ⎥ ⎣ ⎦ (11)

3

⎡ x ⎤ ⎡0 ⎢ y ⎥ ⎢0 ⎢ ⎥=⎢ ⎢ z ⎥ ⎢0 ⎢ ⎥ ⎢ ⎣ w ⎦ ⎣0

1

0

0 0 0 −α 0

ω

0 ⎤ ⎡x⎤ ⎡ 0 ⎤ ⎢ 0 ⎥ 0 ⎥⎥ ⎢⎢ y ⎥⎥ −1 ⎢ ⎥ ( x + bz + aw )3 +P ⎢−K τ 2 ⎥ −ω ⎥ ⎢ z ⎥ ⎥⎢ ⎥ ⎢ ⎥ −α ⎦ ⎣ w⎦ ⎣ 0 ⎦

(14) siendo P-1 la matriz inversa de P dada por la ecuación (13). Teniendo en cuenta las ecuaciones (13), la ecuación (14) en las nuevas coordenadas queda de la forma: ⎡ x ⎤ ⎡ 0 ⎢ y ⎥ ⎢ 0 ⎢ ⎥=⎢ ⎢ z ⎥ ⎢ 0 ⎢ ⎥ ⎢ ⎣ w ⎦ ⎣ 0

1 0 0 0 0 −α 0 ω

0 ⎤ ⎡x⎤ 0 ⎥⎥ ⎢⎢ y ⎥⎥ K + −ω ⎥ ⎢ z ⎥ ω ⎥⎢ ⎥ −α ⎦ ⎣ w⎦

⎡− (ω a + α b ) ⎤ ⎢ ⎥ ⎢ − K pω τ i ⎥ ( x + bz + aw )3 ⎢ ⎥ α ⎢ ⎥ ω ⎣ ⎦

(15) Paso 3. A continuación se simplifica la ecuación (14) tanto como sea posible para reducir la dimensión del sistema a dos. Para ello, y puesto que la no linealidad es cúbica, las coordenadas z y w se aproximan por

⎫ ⎬ 3 3 2 2 w = d1 x + d 2 y + d 3 x y + d 4 xy ⎭ z = c1 x + c2 y + c3 x y + c4 xy 3

3

2

Nótese la acción derivada no aparece en las ecuaciones (21) y no se puede determinar a partir de la bifurcación Bogdanov-Takens. En resumen las condiciones de estabilitdad son: 2α p 2α pf Kp > ⇒ Kp = 2 2 2 2 K α +ω K α +ω (22) 2α pf τi = 2 2 α +ω

2

(16)

donde ci y di (i = 1,2,3,4) se calculan aplicando el teorema de la variedad central. Sin embargo estos coeficientes no son necesarios para analizar la bifurcación, ya que el sistema reducido a la variedad central hasta términos de 5º orden se puede escribir directamente sin tener en cuenta los ci y di, o sea:

⎡ x ⎤ ⎡ 0 1 ⎤ ⎡ x ⎤ ⎡ d1 ⎤ 3 ⎢ y ⎥ = ⎢ 0 0 ⎥ ⎢ y ⎥ + ⎢ d ⎥ x + Ο ( x ⎣ ⎦ ⎣ ⎦⎣ ⎦ ⎣ 2⎦ d1 = −

K

ω

(ωa + α b )

; d2 = −

5

(

(17)

 x = − px 3 ⇒ p

τi

definida positiva. Tomando la derivada respecto del tiempo en la ecuación (18) y sustituyendo las ecuaciones del sistema truncado hasta 5º orden se deduce que 3 3 6 V ( x, y ) = 2 px x + 2 yy = 2 ( p + d 2 ) x y + 2 pd1 x

(19) De las ecuaciones (18) y (19) se deduce que el origen será asintóticamente estable cuando V ( x , y ) < 0 que se cumple si:

)

ω(

τi

(

(20)

=0

)

d1 = − K ω a + α b < 0 ⇒ ω a + α b > 0

ω a + α b = K pω −

(

K α +ω 2

2

)

⇒ Kp >

(

2α p

K α +ω 2

2

G (s) =

)

τ i = KK p

2

(23)

U1 ( s ) U (s)

θ (s) U2 (s)

=

⎛ ⎝

= K p ⎜1 +

1

τis

⎞ ⎠

+τds ⎟

K′

(

s τ s + 2ξτ s + 1 2

2

(24)

)

Suponiendo una entrada senoidal de la forma u1 ( t ) = U 1m sin ωt en el elemento no lineal y que G(s) se asemeja a un filtro paso bajo, la función descriptiva N(U1m) del elemento no lineal viene dada por la ecuación:

(21) en donde se ha considerado la sustitución: K τ = p K. p

+ y 2 = Cte

En la sección anterior la constante τd del regulador PID permanece desconocida, por lo que en esta sección la técnica de la función descriptiva [4] se va a aplicar para investigar la influencia de τd en el comportamiento autooscilante. A la vista de la figura 1 y de las ecuaciones (1) a (9) se deduce que la función transferencia del regulador PID y del conjunto motoramplificador-engranajes y carga viene dadas por las ecuaciones: GPID ( s ) =

siendo p un parámetro dimensional con unidades de V-2s2. De las ecuaciones (17) y (20) se deduce que: 2αω p

x4

4 CALCULO DE LA ACCION DERIVADA CON LA FUNCION DESCRIPTIVA

unidades de V-2s-2 y V ( x, y ) es una función

KK p

)

y por tanto aparece un ciclo límite en el plano de fases del sistema reducido a la variedad central, o sea la condición de auto-oscilación del sistema (f = 1 en la ecuación (23)) se ha determinado. Este comportamiento dinámico es importante puesto que en régimen autooscilante la presencia de cualquier componente armónica en la señal de entrada o en la carga puede dar lugar a comportamiento caótico. En lo que sigue supondremos que p = 1 V-2s2.

en donde los parámetros “a” y “b” viene definidos en la ecuación (13). Una vez que se ha obtenido el sistema reducido a la variedad central, el análisis de estabilidad se lleva a cabo a través del segundo método de Lyapunov. Para ello se escoge una función tentativa de Lyapunov de la forma: 1 4 2 (18) V ( x, y ) = px + y 2 donde p > 0 es un parámetro dimensional con

p + d2 = 0 ⇒ p −

(

en donde el parámetro f ≥ 1 hay que introducirlo para calcular valores admisibles de Kp. De las ecuaciones (17), (20) (21) y (22) tomando f = 1 se deduce que:

)

KK p

)

i

4

3

3 K neU 1m

u2 ( t ) ≈

3

3 K neU 1m

sin ωt ⇒ N (U 1m ) =

4

ω os =

4

(25) De las ecuaciones (24) y (25) la salida G(s) será: θ (t ) =

3

3 K neU 1 m

G ( jω ) sin ( ωt + ϕ1θ )

4 K′

G ( jω ) =

2

ω∆

(

∆ = 4ξ τ ω + 1 − τ ω 2

2

2

(

4ξ τ ω + 1 − τ ω

2

2

)

2

2

2

2

)

2

⎛ 1−τ ω ⎞ = arctg ⎜ ⎟ ⎝ 2ξτω ⎠ 2

2

; ϕ 1θ

2

(26) Teniendo en cuenta las ecuaciones (25) y (26) la salida del regulador PID se puede escribir de la forma: u1 ( t ) =

3

3 K ne K sU 1 m 4

ϕ1 PID = − arctg ⎜

⎝

2

τ iω

f1 min =

⎞ ⎟ ⎠

3

G ( jω ) = U 1 m

4

ϕ1θ + ϕ1 PID + π = 0 ⇒ tg ( ϕ1θ + ϕ1 PID ) = 0 2

U1m =

3 K ne K s G ( jω )

; ω os =

τ i − 2ξτ τ iτ (τ − 2ξτ d )

(28) donde U1m y ωos son la entrada de voltaje en el elemento no lineal y la frecuencia de oscilación respectivamente. De las ecuaciones (28) se deduce:

τ i > 2ξτ ; τ > 2ξτ d ⇒ τ i > 4ξ τ d ⇒ τ d < 2

τi 4ξ

2

⇒ τd <

2 f ξτ 4ξ

2

⇒ τ d = f1

fτ 2ξ

;

(f

1

f −1

τ

f

1 f

⎡ f −1 ⎤ ⎢1 − ⎥ ⎣ f (ω τ ) ⎦ 2

os 1

; f1 max =

1 f

⎡ f −1 ⎤ ⎢1 − ⎥ ⎣ f (ω τ ) ⎦ 2

os 2

2

En esta sección se analizan varios ejemplos para analizar el diseño y comportamiento del servomecanismo controlado con el regulador PID diseñado con el procedimiento indicado en las secciones 3 y 4. El método de simulación utilizado para las ecuaciones (8) y (9) ha sido el de Runge-Kutta de 4º orden, con intervalos de simulación comprendidos entre 0.0001 y 0.001 s. Primero se analizan las condiciones de autooscilación y después se estudia como cambiar los parámetros del regulador PID para alcanzar un punto de consigna deseado. A partir de f = 1 + 10-6 se han obtenido los valores de Kp = 0.0099, τi = 0.0072 s y ωosmin =

(29) y por tanto se ha vuelto a obtener la misma condición para el tiempo de reset τi que se obtuvo en la sección anterior como se puede comprobar de las ecuaciones (21), (22) y (29). Suponiendo que p = 1 y f > 1, el tiempo derivado se puede expresar como: τd <

f (1 − f ⋅ f1 )

1

5 ANALISIS Y DISCUSION DE LOS RESULTADOS OBTENIDOS

τi 4ξ

τ

⇒ ω os min =

(32) Con este procedimiento es posible seleccionar un intervalo de valores para la frecuencia de auto-oscilación. A continuación, seleccionando un valor de f1 tal que f1min < f1 < f1max, es posible determinar la frecuencia de auto-oscilación y el el tiempo derivativo τd del regulador PID. Además si se elige f >> 1, los correspondientes valores de los parámetros Kp y τi del regulador PID pueden llevar al eje de salida del servomecanismo a la posición especificada en el punto de consigna, y por tanto el diseño del regulador PID se ha completado. Hay que resaltar que si el valor de f es muy grande, el tiempo de reset es también grande y por tanto la acción integral es pequeña por lo que un gran error en estado estacionario puede aparecer.

(27) Por otro lado, `puesto que u1(t) es conocida, identificando la amplitud y la fase en las ecuaciones (25)-(27) se deduce que: 3 K ne K sU 1 m

f −1

(31) donde ωosmin es el mínimo valor que puede tomar la frecuencia de auto-oscilación. Hay que tener en cuenta que la frecuencia de autooscilación aparece para f = 1, y por tanto la ecuación (31) no se puede aplicar. Puesto que la constante de tiempo es muy pequeña (ver ecuación (9) y valores de parámetros indicados en la leyenda de la figura 1) se debe escoger f > 1 y muy próxima a la unidad (por ejemplo f = 1.00001 para asegurar muy aproximadamente la condición de auto-oscilación con f.f1 < 1. De los razonamientos anteriores es posible definir un intervalo de frecuencias de auto-oscilación ωos1 > ωosmin y ωos2 > ωos1, y los correspondientes mínimo y máximo valores permitidos para el parámetro f1 o sea:

G ( jω ) sin ( ω t + ϕ1θ + ϕ1 PID + π )

⎛ 1 − τ iτ d ω

1

< 1)

(30) Sustituyendo las ecuaciones (29) y (30) en la ecuación (28) se deduce:

5

frecuencia de oscilación (0.6283 rad/s) está próximo al valor teóricamente calculado con la ecuación (28) (ωos = 0.65 rad/s).

0.2634 rad/s de acuerdo a las condiciones de auto-oscilación dadas por las ecuaciones (22) y (31). A continuación se escogen las frecuencias ωos1 = 0.28 y ωos2 = 3 rad/s con los correspondientes valores f1min = 0.1153 and f1max = 0.9923 calculadas a partir de las ecuaciones (32).Escogiendo una frecuencia de autooscilación deseada de ωos = 0.65 rad/s, se obtiene el valor de f1 = 0.8358 y a partir de la ecuación (30) se obtiene un tiempo derivado de τd = 0.0017. Es importante destacar que el comportamiento auto-oscilante depende de las condiciones iniciales de la señal de control u1(0). Por tanto el valor inicial de las variables de estado se toman todos cero excepto el voltaje x4(0) ≡ u1(0) = U1m = 0.7506 V, que se calcula teniendo en cuenta la ecuación (28), que es la condición de autooscilación obtenida con la función descriptiva.

Figura 3: Comportamiento caótico. Tiempo de simulación t = 600 s, intervalo de simulación T = 0.001 s, Ka = 100, ξ = 0.95, τ = 3.7969.10-3 s, Km = 4.3373 N.m/A, p = 1, f = 1 + 10-5, ωosmin = 0.8328 rad/s, ωos1 = 0.84, ωos2 = 3, f1min = 0.017, f1max = 0.9229, ωos = 1 rad/s, f1 = 0.3064, Kp = 0.0099, τi = 0.0072 s, τd = 6.1224.10-4 s, U1m = 1.1547, entrada de referencia en tr = 100 s, θr1 =1, θr2 = 3.456, ωr = 0.987 ras/s y Td2 = 0. Condiciones iniciales xo = [0 0 0 1.1547]. Cambio valores regulador PID en tcon = 400 s con f = 200, Kpc = 1.9718, τic = 1.4428 s y τdc = τd = 6.1224.10-4 s, punto de consigna θsp = -0.5 rad, región de captura: rcx = 0.8 V, rcy = 0.2 V/s. La figura 3 muestra un comportamiento dinámico muy interesante cuando se varía armónicamente el punto de consigna a partir del instante t = 100 s (ver ecuaciones 8 y 9). Como en el caso de la figura 2, se eligen valores de Kp, τi y τd para obtener comportamiento autooscilante con ωos = 1 rad/s (el valor obtenido en la simulación es 0.9844 rad/s). En el instante t = 100 s el punto de consigna se varía de cero a θ ( t ) = 1 + 3.456 sin ( 0.987t ) , lo que da lugar a

Figura 2: Tiempo de simulación t = 200 s, , Ka = 100, ξ = 0.95, τ = 3.7969.10-3 s, Km = 4.3373 N.m/A, p = 1, f = 1 + 10-6, ωosmin = 0.2634 rad/s, ωos1 = 0.28, ωos2 = 3, f1min = 0.1153, f1max = 0.9923, ωos = 0.65 rad/s, f1 = 0.8358, Kp = 0.0099, τi = 0.0072 s, τd = 0.0017 s, U1m = 0.7506, θr1 = θr2 = 0 , Td2 = 0. Condiciones iniciales xo = [0 0 0 0.7506]. El control para estabilidad se aplica en tcon = 50 s with f = 200, Kpc = 1.9718, τic = 1.4428 s and τdc = τd = 0.0017 s, punto de consigna θsp = 0.8 rad.

r

la aparición de oscilaciones caóticas como se muestra en la figura 3 a). En la figura 3 b) se muestran dos valores de la variable x1(t) obtenidos con dos condiciones iniciales que difieren en 10-10. Se observa que para t ≈ 250 s las oscilaciones son completamente diferentes con fuerte dependencia sensible, que es un indicador de comportamiento caótico. [2], [3], [13]. Puesto que el punto de consigna está dentro del atractor extraño, que aparece como consecuencia del comportamiento caótico, siempre existirá una órbita caótica que pase muy cerca del punto de consigna. Por consiguiente, cambiando el valor del parámetro f (ecuaciones (22) y (30)) cuando una órbita

Los resultados de la simulación y los valores de los parámetros utilizados se muestran en la figura 2. En las figuras 2 a) y 2 b) se muestran ciclos límites durante el tiempo 0 ≤ t ≤ tcon ya que los parámetros f ≈ 1 y f1 = 0.8358 dan valores de las constantes del regulador PID que dan comportamiento auto-oscilante. En el instante tcon = 50 se cambia f = 200 y se obtienen nuevos valores de Kp y τi tal que el punto de equilibrio es asintóticamente estable. La figura c) muestra la respuesta temporal y en la d) se ha calculado la densidad espectral de potencia de la señal u1(t) para comprobar que la

6

cruce una pequeña zona de captura (rcx, rcy) definida alrededor del punto de consigna (ver figura 3 c)) será posible posicionar el servomecanismo en un punto de consigna predeterminado, con valores muy pequeños de la señal de control u1(t). Puesto que el momento en que esto suceda es a priori desconocido, la región de captura debe de elegirse para evitar largos tiempos de espera antes de alcanzar el punto de consigna. Los valores de los parámetros para los diferente comportamientos dinámicos se indican en la leyenda de la figura 3.

f para obtener estabilidad asintótica para el punto de consigna elegido que en este caso es θsp = 0.5 rad. En un instante arbitrario de tp = 10 s se aplica un par de perturbación externo de la forma: T ( t ) = 10.987 sin ( 3.456t ) N y el sistema d

salta a un ciclo límite y se mantiene en él hasta que en el instante t = tcon = 30 s el valor del parámetro f se aumenta hasta f = 150. Suponiendo la presencia de ruido definido por las ecuaciones (32) con fna = 0.2 se comprueba que se alcanza el punto de consigna, aún en presencia de ruido.

6 CONCLUSIONES En este trabajo se ha estudiado el comportamiento dinámico de un servomecanismo de posición con una no linealidad cúbica controlado por un regulador PID. Se pretende estabilizar un cuerpo conectado al eje de salida del servomecanismo que está sometido a un par de perturbación externo. Debido a la presencia de la no linealidad, se demuestra que en el sistema realimentado con el regulador PID, aparece una bifurcación del tipo Bogdanov-Takens, con una matriz de la parte lineal del sistema con dos autovalores cero. El análisis de la bifurcación, calculando todos los autovalores de la parte lineal, junto con el teorema de la variedad central, permite deducir un sistema simplificado de orden dos, que posee las mismas propiedades dinámicas del sistema original. Este sistema se usa para analizar la estabilidad a partir del segundo método de Lyapunov, lo que permite definir la constante de acción proporcional e integral del regulador PID de forma que el punto de equilibrio sea asintóticamente estable. El tiempo derivado τd junto a la frecuencia de auto-oscilación se han determinado a través de la función descriptiva de la no linealidad cúbica. Además se ha comprobado que variaciones armónicas en la entrada de referencia pueden producir oscilaciones caóticas, las cuales se pueden usar junto con el regulador PID para obtener señales de control muy pequeñas El correcto funcionamiento del sistema realimentado se ha comprobado suponiendo una perturbación armónica en el eje de salida del servomecanismo que está afectada por ruido aleatorio. El análisis global del sistema da una visión unificada en donde el control PID, la bifurcación de Bogdanov-Takens, la función descriptiva, el segundo método de Lyapunov y el comportamiento caótico muestran un procedimiento de cómo fenómenos complejos que no son fáciles de interpretar pueden tratarse de forma unificada.

Figura 4:: Tiempo de simulación t = 40 s, intervalo de simulación T = 0.004 s, Ka = 10, ξ = 0.95, τ = 3.7969.10-3 s, Km = 4.3373 N.m/A, p = 1, f = 10, Kp = 0.9895, τi = 0.0721 s, τd = 0.004 s, θr1 = 1, θr2 = 0, par de perturbación aplicado en tp = 10 s, Td2 = 10.987 N y ωp =3.456 rad/s. Condiciones iniciales xo = [0.3183 0 0 0]. Cambio de valores del regulador PID en tcon = 30 s con f = 150, Kpc = 14.7881, τic = 1.0821 s y τdc = τd = 0.004 s. Punto de consigna θsp = 0.5 rad factor de ruido fna = 0.2. Para comprobar la validez del diseño, se supone que en el sistema puede estar afectado por ruido aleatorio n(t) tal como se indica en la figura 1 y por tanto las variables de estado toman los valores x1 ( t ) = x1 ( t ) + f na [ X − 0.5]

x2 ( t ) = x2 ( t ) + f na [ X − 0.5] x3 ( t ) = x3 ( t ) + f na [ X − 0.5]

(33)

x4 ( t ) = x4 ( t ) + f na [ X − 0.5]

donde X es una variable aleatoria que está uniformemente distribuida entre 0 1 y fna > 0 es un factor de amplificación para obtener una amplitud de ruido distribuida uniformemente entre –fna/2 y fna/2. Como se muestra en la figura 4, primero reeligen valores de Kp, τi, τd y 7

La metodología presentada en este trabajo se puede extender a sistemas con no linealidades más complejas y con estrategias de control más sofisticadas que pueden dar lugar a comportamientos dinámicos más interesantes. Referencias

[1] Bogdanov R., Versal deformations of a singular point on the plane in the case of zero eigenvalues, Selecta Math. Soviet.1 (1981) 389–421. [2] Chiba Hayato, Periodic orbits and chaos in fast–slow systems with Bogdanov Takens type fold points, J. Differential Equations. 250 (2011) 112-160. [3] Guckenheimer J., Holmes P., Nonlinear Oscillations, Dynamical Systems and Bifurcations of Vector Fields. Springer Verlag, New York, 1983. [4] Khalil K., Nonlinear Systems, second ed., Prentice Hall, New Jersey, 1996. [5] Kuznetsov Y. A., Elements of Applied Bifurcation Theory, third ed., Springer Verlag, New York, 1998. [6] Péréz. Polo M, Albertos P.,.Berná Galiano J. A., Tuning of a PID controller gyro using the bifurcation theory, Syst. Control Lett. 57 (2008) 10-17. [7] Pérez Polo M. F. Pérez. Molina M., Gil Chica J., Stabilization and positioning control of a rolling disk by using the Bogdanov-Takens bifurcation, Physica D. 241 (2012) 1450-1469. [8] Pérez-Molina M., Pérez-Polo M. F., Steady-state, self-oscillating and chaotic behavior of a PID controlled nonlinear servomechanism by using the BogdanovTakens and Andronov-Poincaré_Hopf bifurcations, Commu. Nonlinear Sci. Numer. Simulat., 19 (2014), 3694-3717. [9] Starzhinskii V. M., Applied Methods in the Theory of Nonlinear oscillations. Mir Publishers, Moscow, 1980 [10] Takens F., Forced oscillations and bifurcations, Comm. Math. Inst. Rijksuniv. Utrecht. 3 (1974) 1–59. [11] Verduzco F., Alvarez J., Hopf bifurcation control: A new approach, Syst. Control Lett. 55 (2006) 437-451. [12] Verduzco F., Carrillo F. A., TakensBogdanov Bifurcation Analysis in Indirect Field-oriented Control, AIP Conference Proceedings. 1368 (2011) 101-104 [13] Wiggins S., Introduction to Applied Nonlinear Dynamical Systems and Chaos, second ed., Springer Verlag, New York, 1990.

8