Ecuaciones en Diferencias

Ecuaciones en diferencias lineales Félix eli x Monas Mo naster terioio-Hue Huelin lin 11 de abril de 2010

´ Indice 1. Introducci Introducci´ on o ń

2

2. Estudio de la estabilidad

9

2.1. Teorema de Routh-Hurwitz . . . . . . . . . . . . . . . . . . . 12 3. La transfor transformada mada

Z

16

4. Resol Resoluci uci´ on o ń de la ecuaci´ on on en diferencias diferenc ias utilizando utiliz ando el m´ etodo etod o de la descomposici´ on on o expansi´ on o n en frac fracci cion ones es parc parcia iale les. s. 20

4.1. 4.1. Caso Caso de polos reales reales simples simples,, con al menos un cero en el or´ or´ıgen 22 4.2. 4.2. So Solu luci ci´ón o n con con po pollos comp comple lejo joss simp simple less . . . . . . . . . . . . . . 25 4.3. Soluci´ ción con po pollos dob oblles . . . . . . . . . . . . . . . . . . . . . 29 5. Resol Resoluci uci´ on o ń de la ecuaci´ on en diferencias utilizando el teorema on de los residuos 32

1

1.

Introd trodu ucci´ cci´ on on

En estos breves apuntes voy a limitarme al an´ alisis de ecuaciones en diferalisis encias lineales de parámetros ametros constantes. Supondré tambi´ en en que todos los parámetros ametros son so n conocidos, cono cidos, es decir que no haré ninguna ningun a menci´ on on a las técniecn icas que se hayan utilizado para obtenerlos. En algunos casos estas ecuaciones pueden provenir de una an´ alisis alisis estad´ estad´ıstico de series temporales, en otros de modelos basados en leyes leyes cient cient´ıficas y en otros de la discretizaci´ discretizaci´ on de ecuaciones ecuaciones diferenciales diferenciales continuas. continuas. En este breve estudio el or´ or´ıgen de las ecuaciones es indiferente ya que voy a centrar el an´ alisis a lisis en el estudio de propiedades propiedades generales de las ecuaciones ecuaciones y no en el significado significado f´ısico de las variables y parámetros. ametros. No obstante o bstante supondr´ sup ondré principalmente que estas ecuaciones representan variaciones en un contexto temporal. 1. Entradas y salidas Tambi´ ambién en voy voy a limitar limitar el estudio al caso de una unica uńica variable, que lla ll amaré y, dependiente impl´ impl´ıcitamente de un ´ındice de variaci´ on k y que en general será funci´ on lineal de una, ninguna o varias variables. on En términos erminos de sistemas a la variable y la llam ll amar ar´é salida del sistema, y a las restantes, entradas. Por ejemplo la ecuaci´ on on yk + 2 yk−1 = 3uk−1 + 5 vk−1

representa un sistema cuya entradas son u y v y cuya salida es y. 2. Principio de superposici´ on on Podemos observar que la ecuaci´ on es lineal, lo que significa que para su on resolución on puede utilizarse el principio de superposici´ on on. Esto nos permite resolver la ecuaci´ on anterior resolviendo independientemente on dos ecuaciones, ecuaciones, ykv + 2ykv−1 = 5vk−1 yku + 2yku−1 = 3uk−1 ya que la linealidad del sistema garantiza que la salida será la suma de las salidas parciales debidas a cada una de las entradas por separado, yk = ykv + yku

2

1.

Introd trodu ucci´ cci´ on on

En estos breves apuntes voy a limitarme al an´ alisis de ecuaciones en diferalisis encias lineales de parámetros ametros constantes. Supondré tambi´ en en que todos los parámetros ametros son so n conocidos, cono cidos, es decir que no haré ninguna ningun a menci´ on on a las técniecn icas que se hayan utilizado para obtenerlos. En algunos casos estas ecuaciones pueden provenir de una an´ alisis alisis estad´ estad´ıstico de series temporales, en otros de modelos basados en leyes leyes cient cient´ıficas y en otros de la discretizaci´ discretizaci´ on de ecuaciones ecuaciones diferenciales diferenciales continuas. continuas. En este breve estudio el or´ or´ıgen de las ecuaciones es indiferente ya que voy a centrar el an´ alisis a lisis en el estudio de propiedades propiedades generales de las ecuaciones ecuaciones y no en el significado significado f´ısico de las variables y parámetros. ametros. No obstante o bstante supondr´ sup ondré principalmente que estas ecuaciones representan variaciones en un contexto temporal. 1. Entradas y salidas Tambi´ ambién en voy voy a limitar limitar el estudio al caso de una unica uńica variable, que lla ll amaré y, dependiente impl´ impl´ıcitamente de un ´ındice de variaci´ on k y que en general será funci´ on lineal de una, ninguna o varias variables. on En términos erminos de sistemas a la variable y la llam ll amar ar´é salida del sistema, y a las restantes, entradas. Por ejemplo la ecuaci´ on on yk + 2 yk−1 = 3uk−1 + 5 vk−1

representa un sistema cuya entradas son u y v y cuya salida es y. 2. Principio de superposici´ on on Podemos observar que la ecuaci´ on es lineal, lo que significa que para su on resolución on puede utilizarse el principio de superposici´ on on. Esto nos permite resolver la ecuaci´ on anterior resolviendo independientemente on dos ecuaciones, ecuaciones, ykv + 2ykv−1 = 5vk−1 yku + 2yku−1 = 3uk−1 ya que la linealidad del sistema garantiza que la salida será la suma de las salidas parciales debidas a cada una de las entradas por separado, yk = ykv + yku

2

En consecuencia, la utilizaci´ on del principio de superposici´ on on on nos permite reducir el estudio a sistemas de una unica u ´ nica entrada y una unica uńica salida, sin ninguna pérdida erdida de generalidad. 3. Principio de causalidad La terminolog terminolog´´ıa de entrada-sali entrada-salida da tiene un sentido sentido especial cuando se está en un contexto temporal, y sobre todo si se admite el principio de causalidad. Bajo este supuesto las salidas en un instante de tiempo no pueden depender de entradas en instantes de tiempo futuros, y es por esta razón on que tiene sentido decir, en el ejemplo anterior, que y es la salida y u y v son entradas. Un ecuaci´ on on de la forma yk + 2 yk−1 = 3uk+1

nos obligar´ obligar´ıa a decir que u es la salida e y la entrada, si se admite el principio de causalidad. Esto significa que bajo este principo no hay ninguna ambig¨ uedad uedad en el concepto de entrada y salida, excepto cuando el orden sea el mismo para las dos variables. Por ejemplo, yk + 2 yk−1 = 3uk

no nos permite saber si u o y son las entradas o salidas. En términos erminos de sistemas este ultimo u´ltimo hecho tiene un significado f´ısico importante: la salida responder´ a a la entrada instant´ aneamente. aneamente. En estos apuntes utilizaré la notaci´ on on y para representar la salida, y representar la entrada, entrada, y admitir´ admitiré que pueda haber respuestas respuestas u para representar intantáneas, aneas, pero p ero a la vez supondr´ sup ondré que siempre se cumple el principio de causalidad. Cuando el orden de la entrada sea inferior al de salida diré que se cumple c umple un principio de d e causalidad ca usalidad estricta. est ricta. 4. Or´ Or´ıgen de tiempos y condiciones iniciales Esta imposición on presenta, no obstante, un problema relacionado con el or´ est e es or´ıgen de tiempo tiem poss y con las condiciones iniciales. Pero éste un problema menor si se toma en consideraci´ on que siempre es posible on realizar una traslaci´ on en el tiempo. La ecuación on on yk + 2 yk−1 = 3uk−1

3

puede convertirse, sin pérdida de generalidad, en la ecuaci´ on yk+1 + 2 yk = 3uk

quedando el or´ıgen de tiempos situado en k = 0, y las condiciones inciales en y0 , donde en general y0 = 0.



En la ecuaci´ on regresiva la condici´ on inicial es y−1 que en general es distinto de cero, por lo tanto a la hora de interpretar las condiciones iniciales deberá prestarse atenci´ o n a la forma en que se represente la ecuación en diferencias, ya que en una u otra se est´ a considerando un or´ıgen de tiempos distinto. O dicho de otra forma, k = 0 no define de manera general el or´ıgen de tiempos, salvo que se expresen las ecuaciones en su forma causal. En la forma causal se presupone que los valores de las entradas y salidas anteriores al or´ıgen de tiempos son nulas. Esto da sentido al concepto de or´ıgen de tiempos, y como he mostrado pueden evitarse errores de interpretaci´ on si se hace que k = 0 represente este or´ıgen de tiempos. No obstante, matem´ aticamente no aparece ning´ un problema con una u otra representaci´ on. En lo que sigue utilizaré indistintamente una u otra representaci´ on, con la precaución mencionada de saber cu´ al es el or´ıgen de tiempos. Desde un punto de vista matem´ atico el ´ındice de variación k no tiene por qué representar el tiempo, sino cualquier otra cosa. No obstante en este estudio supondré que k 0, 1, 2 . . . , que se corresponder´ a con el tiempo f´ısico t = kT donde T es un n´ umero fijo y conocido expresado en segundos, a˜ nos o cualquier otra unidad de tiempo, que llamaré periodo de muestreo. Realmente no es necesario hacer este supuesto temporal, pero resultar´ a conveniente m´ as adelante cuando quieran interpretarse los resultados en un contexto temporal. Pero como digo, no es necesario dar al ´ındice k ning´ un significado f´ısico especial; es suficiente con considerarlo un ´ındice de variaci´ on abstracto.

∈{

}

5. Ecuación homogénea Llamaré ecuaci´ on homogénea de una ecuación en diferencias o de un sistema, a aquella cuya entrada sea nula. Por ejemplo, la ecuación homogénea del sistema yk + 2 yk−1 = 3uk−1

4

es yk + 2yk−1 = 0

La solución de la ecuaci´ on homogénea depende de las condiciones iniciales. Si éstas son nulas, entonces yk = 0. Escribiendo la anterior ecuaci´ on en su forma causal, yk+1 + 2yk = 0

podemos ver que y1 = ción general

−2y , pudiéndose obtener por inducción la soluy = (−2) y , k = 0, 1, 2, . . . k

0

k

0

Para la resolución completa de la ecuación yk + 2 yk−1 = 3uk−1

es necesario conocer la funci´ on uk , k. Puede demostrarse que si esta función es conocida y se obtiene una soluci´ on particular para yk , la solución completa se obtiene como suma de la soluci´ on particular y de la solución general de la ecuaci´ on homogénea. La soluci´ on particular se obtiene suponiendo condiciones iniciales nulas, es decir y0 = 0 en el ejemplo anterior.

∀

Veremos en las siguientes secciones c´ omo obtener la soluci´ on particular. 6. Representaci´ on en el espacio de estados Hasta aqu´ı se ha representado un sistema como una ecuaci´ on en diferencias a través de sus entradas y salidas. Por ejemplo, el siguiente es un sistema de orden dos con una entrada u y una salida y , yk + 2 yk−1

− 5y −

= 3uk−2

k 2

que en la forma causal puede escribirse como yk+2 + 2 yk+1

− 5y

k

= 3uk

Para la resoluci´ on de este sistema es necesario conocer, tanto la entrada, como las condiciones iniciales de la salida, y0 e y1 .

5

Vamos a ver en este apartado que realizando una transformaci´ on adecuada, toda ecuaci´ on en diferencias de orden n es equivalente a n ecuaciones en diferencia de orden 1, que llamaremos ecuaci´ on de estados, junto con una ecuaci´ on adicional que llamaremos ecuaci´ on de salida. Llamemos

x1,k = yk x2,k = yk+1

Vemos en primer lugar que x2,k = x1,k+1

y en segundo lugar que yk+2 = x2,k+1

Por lo tanto x2,k+1

x1,k+1 = x2,k = 2x2,k + 5x1,k + 3uk

−

obteniendo as´ı dos ecuaciones en diferencias de primer orden, a través del vector de estados xk =

  x1,k x2,k

Estas dos ecuaciones pueden escribirse en forma matricial



x1,k+1 x2,k+1

 0 =

5

1 2

−

   0  x1,k x2,k

+

3

uk

En forma más compacta obtendr´ıamos la ecuaci´ on matricial llamada ecuación de estados, xk+1 = Axk + Bu k donde A y B son matrices de dimensiones [n n] y [n 1] respectivamente, siendo n el orden del sistema, y xk el vector de estados, de dimensión [n 1].

×

×

×

A esta representaci´ on se le llama representación en el espacio de estados o también representaci´ on interna, mientras que a la anterior se la llama representaci´ on de entrada-salida.

6

La ecuaci´ on de estados no define completamente un sistema de entradasalida. Es necesario especificar la relaci´ on del vector de estados con la salida. Esto se hace a través de la ecuaci´ on de salida yk

    = 1 0 x1,k x2,k

o en forma compacta yk = Cxk

donde C es una matriz de dimensión [1

× n].

Podemos observar que las condiciones iniciales en esta representaci´ on vienen dadas por el vector x0 , ya que x0 =

  y0 y1

En estos breves apuntes no voy a incluir có mo se obtiene la representaci´ o n en el espacio de estados de una ecuación en diferencias de orden n cualquiera, aunque se ha visto una forma de hacerlo para un caso particular. Concretamente para el caso en que la expresión de la entrada dependa exclusivamente de uk . Por el momento es suficiente con saber que cualquier ecuaci´ on en diferencias lineal de orden n puede expresarse como una ecuaci´ on de estados de la forma xk+1 = Axk + Bu k , junto con una ecuaci´ on de salida, que no es una ecuaci´ on en diferencias sino una relaci´ o n directa entre la salida y el vector de estados. Si la ecuación en diferencias es lineal y de parámteros constantes, las matrices A y B serán constantes. El interés que tiene esta representaci´ on es muy variado, aunque desde el primer momento puede observarse que para la resolución de ecuaciones en diferencias de cualquier orden es suficiente con saber resolver ecuaciones en diferencias matriciales de primer orden. Por ejemplo, la solución de la ecuación homogénea es casi trivial, ya que si xk+1 = Axk

entonces xk = Ak x0 ,

k = 0, 1, 2 . . .

y la solución general de la ecuaci´ on homogénea de la salida ser´ a yk = CAk x0 ,

7

k = 0, 1, 2 . . .

Otro aspecto resaltable de esta representaci´ o n es que el vector de es´ tados tiene una interpretaci´ on f´ısica interesante. Este expresa toda la información necesaria para describir el comportamiento de un sistema, además de que muestra con claridad que un sistema de orden n puede entenderse como un objeto geométrico de dimensi´ on n, en el que cada estado xi representa una coordenada. En el caso bidimensional puede utilizarse un plano de coordenadas x = (x1 , x2 ), llamado plano de fase, y dibujarse la evolución del vector o punto xk variando el ´ındice k, desde k = 0, es decir desde x0 , hasta k = . Hablaré de esto m´ as adelante.

∞

7. Comentario sobre los procesos estocásticos Un sistema cuyas variables sean variables aleatorias se le llama proceso o sistema estoc´ astico. Por ejemplo un sistema determin´ıstico que tenga una entrada aleatoria se transforma en un sistema estocástico. Por ejemplo, si el sistema determin´ıstico yk + 2 yk−1 = 3uk−1

tiene una entrada adicional vk aleatoria, como por ejemplo un ruido blanco o una entrada que siga una distribución gaussiana, que influya en el sistema de la siguiente forma, yk + 2 yk−1 = 3uk−1 + 5 vk−1

transformará el sistema determin´ıstico en un sistema o proceso estocástico. A partir de este momento el estudio del sistema requiere utilizar la teor´ıa de los procesos estoc´ asticos. Por ejemplo, si la entrada vk tiene media cero o constante, es decir si E [vk ] = c, entonces puede trabajarse con valores medios (E representa realmente la esperanza matem´ atica) como si fuese el sistema determin´ıstico, E [yk+1] + 2 E [yk ] = 3uk + 5 c

La entrada uk sigue siendo determin´ıstica, pero la salida es ahora una variable aleatoria. Puesto que en este sistema no se encuentra toda la informaci´ o n estad´ıstica será preciso incorporarla de alguna otra forma, como por ejemplo utilizando expresiones de la covarianza. No voy a hablar en estos apuntes de esta clase de sistemas. 8

2.

Estudio de la estabilidad

En estos apuntes voy a hacer varios estudios. Por un lado la obtención de la solución completa de una ecuaci´ o n en diferencias, y por otro el an´ alisis de esta solución. Un estudio fundamental ser´ a saber si el sistema es estable. Realmente el concepto de estabilidad no se aplica a los sistemas (si se refiere a éstos se suele hablar de estabilidad estructural) sino a sus puntos de equilibrio (que podr´ıa llamarse estabilidad funcional), pero los sistemas lineales tienen un u ´ nico punto de equilibrio por lo que habitualmente se habla de estabilidad del sistema abusando de lenguaje. También suele haber una confusi´ on en el concepto de sistema, ya que cuando se dice que un sistema lineal tiene un único punto de equilibrio se sobreentiende que se est´ a refiriendo a un sistema cuyas entradas son nulas, es decir a la ecuación homogénea del sistema. La entradas pueden ser muy diversas, y pueden desestabilizar un sistema lineal o convertirlo en un sistema no lineal. Para el estudio de la estabilidad es suficiente con estudiar la matriz de estados A, o equivalentemente la ecuaci´ on caracter´ıstica. Dado un sistema expresado en el espacio de estados, xk+1 = Axk , se definen los puntos de equilibrio como aquellos que cumplen que xk+1 = xk = x∗ . Es decir x∗ = Ax∗ , o también

− A)x∗ = 0 Puesto que en general, el determinante de (I − A) es distinto de cero, la uńica ∗ (I

solución posible es que x = 0, es decir que los sistemas lineales tienen un u ´ nico punto de equilibrio, que es precisamente el nulo. Solamente en el caso en que la matriz A tenga alg´ un autovalor unidad, puede pensarse que un sistema lineal tenga infinitos puntos de equilibrio. Esta situaci´ on se analizar´ a como un caso particular. Con este simple análisis hemos calculado el punto de equilibrio, pero no se ha estudiado su estabilidad. Se define la ecuación caracter´ıstica del sistema como det(λI

− A) = 0 9

es decir, que las soluciones de la ecuaci´ on caracter´ıstica son los autovalores de la matriz A. Puede demostrarse que el sistema es estable si los m´ odulos de sus autovalores son menores que la unidad. Pero antes de avanzar deberemos saber qu´ e se entiende por estabilidad de un punto de equilibrio. Hay una buena colecci´ on de definiciones de estabilidad, si bien casi todas responden a la idea intuitiva de Lyapunov: si un sistema está en un punto de equilibrio y éste se perturba de alguna forma (con una peque˜ na perturbación), el sistema tender´ a a regresar a él si el punto de equilibrio es estable, o tender´ a a alejarse de él si se trata de un punto de equilibrio inestable. En el contexto de la teor´ıa de sistemas hay otra idea intuitiva de estabilidad: si la entrada de un sistema est´ a acotada y observamos que su salida también está acotada, entonces el sistema ser´ a estable. Ambas ideas sirven como definici´ on de estabilidad y puede comprobarse que para los sistemas lineales son equivalentes. Puesto que los sistemas lineales tienen un uńico punto de equilibrio, cualquier punto del espacio puede considerarse una perturbación del punto de equilibrio. Por lo que el estudio de la estabilidad se reduce a averiguar si, para cualquier condici´ on inicial, el estado se acercará al punto de equilibrio o se alejará de él. En los sistemas no lineales puede haber muchos puntos de equilibrio, por lo que deber´ a estudiarse con cuidado la clase de perturbaciones que puede realizarse, y es por esto que es necesario calcular las regiones de atracción o de repulsión de cada uno de los puntos de equilibrio. Pero en los sistemas lineales, todo el espacio o es de repulsión o es de atracción. En el ejemplo del anterior apartado hab´ıamos visto que la matriz de estados de la ecuación en diferencias yk+2 + 2yk+1

es A=

− 5y

0 5

por lo que la ecuaci´ on caracter´ıstica es det(λI

k

1 2

−

= 3uk



− A) = λ(λ + 2) − 5 = 0 10

que podemos escribir en la forma λ2 + 2 λ

−5 =0

Podemos observar que se trata de un polinomio cuyo orden es el mismo que el orden del sistema, y cuyos coeficientes son los coeficientes del sistema homogéneo de la ecuación en diferencias. Puede demostrarse que esto siempre es as´ı.

− ±√

En este ejemplo las ra´ıces del polinomio caracter´ıstico son 1 6, por lo que como ambas ra´ıces son en m´ odulo mayores que la unidad, el sistema es inestable. Bastar´ıa con que una de ellas lo hubiese sido para que el sistema fuese inestable. Puede comprobarse tambi´ en que la soluci´ on homogénea xk = Ak x0 tiende a infinito (en norma 1 ) cuando k , para cualquier valor de las condiciones iniciales x0 , lo que indica que el estado se ha alejado del punto de equilibrio. Sin embargo, si el sistema hubiese sido estable se cumplir´ıa que limk→+∞ xk = 0, es decir que el estado tender´ıa al punto de equilibrio, para cualquier condici´ on inicial. Por otro lado, puesto que la salida del sistema se relaciona con el vector de estados a través de una matriz constante yk = Cxk , se deber´ a cumplir que si el sistema es inestable lim k→+∞ yk = , y si el sistema es estable limk→+∞ yk = 0.

→∞

|| ||

±∞

En este estudio de la estabilidad no se ha visto el efecto de las entradas. Como ya he se˜ nalado su efecto dependerá de qué clase de entradas haya. ´ Estas podr´ an desestabilizar un sistema estable o estabilizar uno inestable. Por lo tanto la estabilidad de un sistema puede considerarse una propiedad intr´ınseca al sistema. No obstante puede garantizarse que si las entradas están acotadas, como es el caso de entradas constantes, por ejemplo, éstas no afectar´ an a la estabilidad o inestabilidad del sistema. Antes de hacer otros estudios conviene prestarle atenci´ o n al hecho de que para saber si un sistema es estable o no, es necesario calcular las ra´ıces de un polinomio de orden n. Se sabe por el teorema fundamental del a´lgebra que la ecuación caracter´ıstica tendrá n ra´ıces, y que éstas ser´ a n reales o complejas, con un orden de multiplicidad cualquiera; tambi´ en se sabe que si 1 Se

cumple que xk A k x0 , y si los autovalores de A son mayores que la unidad en módulo, A k tiende a infinito cuando k .

|| ||

|| ||≤|| || || ||

→∞

11

son complejas siempre aparecer´ an en pares conjugados. Pero tambi´ en se sabe que en general, no es posible resolver el polinomio por radicales cuando n 5, lo que plantea un problema numérico que puede tener importancia pr´ actica. Se conocen muchos métodos numéricos que dan buenas aproximaciones a las ra´ıces de una ecuaci´ on, pero debido a que no son soluciones exactas pueden distorsionar el an´ alisis de la estabilidad.

≥

Vamos a ver a continuación un método que permite saber si un sistema es estable o no sin necesidad de calcular expl´ıcitamente las ra´ıces de la ecuaci´ on caracter´ıstica.

2.1.

Teorema de Routh-Hurwitz

Hay diversas formas de saber si las ra´ıces de un polinomio caen dentro del c´ırculo unidad del plano complejo (es decir, con m´ odulo menor que la unidad). El plano complejo se define como se hace habitualmente, es decir siendo el eje de abcisas la parte real (σ ) y el eje de ordenadas la parte imaginaria (ν ) de un n´ umero complejo w = σ + iν . El módulo ser´ a w 2 = σ2 + ν 2 , que representa una circunferencia de radio w en el plano complejo.

| |

| |

Un teorema directamente aplicable a este estudio de estabilidad es el teorema de Jury, y otro es el teorema de Routh-Hurwitz. Voy a explicar el segundo método que aunque fue originalmente concebido para el estudio de la estabilidad de sistemas continuos representados mediante ecuaciones diferenciales, puede aplicarse al caso de sistemas discretos representados mediante ecuaciones en diferencias. Dada la ecuación caracter´ıstica (que escribiré utilizando la variable compleja z en vez de λ como hice anteriormente), en la forma siguiente: a0 z n + a1 z n−1 + a2 z n−2 +

expresada de tal forma que a0 > 0.

··· + a

n

=0

En primer lugar se realiza la transformación bilineal z=

w+1 w 1

−

obteniendo otro polinomio que tendr´ a la forma b0 wn + b1 wn−1 + b2 wn−2 +

12

··· + b

n

=0

donde de nuevo la escribimos de tal manera que b0 > 0. Puede demostrarse que esta nueva ecuación caracter´ıstica expresada en términos de w tiene ra´ıces cuya parte real es negativa si y solamente si la ecuaci´ on expresada en términos de z tiene ra´ıces dentro del c´ırculo unidad. Esto puede demostrarse analizando la transformaci´ on bilineal entre z y w. Por ejemplo, para el sistema que se ha estudiado hasta ahora, z 2 + 2z

−5=0

se obtiene +1 2 +1 (w ) + 2( w ) w −1 w −1

−5 =0

y de aqu´ı, (w + 1)2 + 2(w + 1)(w

2

− 1) − 5(w − 1)

=0

resultando w2

− 6w + 3 = 0

Puede considerarse que esta nueva ecuaci´ on se corresponde con la ecuaci´ on caracter´ıstica de un sistema continuo, para el cual la estabilidad queda garantizada si la parte real de las ra´ıces es menor que cero. Queda claro en este ejemplo que el sistema continuo equivalente ser´ a inestable ya que las ra´ıces son 3 6, que son positivas. En consecuencia lo será el sistema discreto del que deriva.

±√

Para la aplicaci´ on del teorema de Routh-Hurwitz es necesario construir en primer lugar una tabla de (n + 1) filas, cuyas dos primeras filas se forman con los coeficientes de la ecuación caracter´ıstica. En la primera fila se colocan los coeficientes de orden par si n es par o los de orden impar si n es impar. La segunda fila se forma con los restantes coeficientes. Las siguientes filas se calculan como se ver´ a en los ejemplos. Posteriormente se estudian los cambios de signo de su primera columna. Si no hay cambios de signo, el teorema garantiza que no hay ninguna ra´ız inestable; además cada cambio de signo indica la existencia de una ra´ız inestable. Veámoslo en primer lugar en el ejemplo anterior, w2

− 6w + 3 = 0 13

La tabla de Routh-Hurwitz es: w2 w1 w0

1 3 6 0

−

c

donde c se calcula de la siguiente forma: c=

( 6)(3) (1)(0) ( 6)

−

− −

=3

La columna de la izquierda no tiene ning´ un significado especial. Sirve tan sólo para indicar el orden de los polinomios que se van construyendo. Vemos que hay dos cambios de signo en la primera columna (de n´ umeros), entre 1 y 6 y entre 6 y 3, por lo que el teorema de Routh-Hurwitz garantiza que la ecuación caracter´ıstica tiene dos ra´ıces inestables, es decir con parte real positiva. En consecuencia el sistema dado por

−

−

z 2 + 2z

−5=0

tendrá dos ra´ıces fuera del c´ırculo unidad. Veamos otro ejemplo de un sistema de orden 5, sobre el que ya se ha realizado la transformaci´ on bilineal, obteniendo el polinomio w5 + 2w4 + 24w3 + 50w2

− 25w − 52 = 0

La tabla de Routh-Hurwitz es: w5 1 24 w4 2 50 w3 c0 c1 w2 d0 d1 w1 e w0 f

−25 −52 0

Las restantes filas se calculan como se hizo en el anterior ejemplo. As´ı para la fila w3 será (2)(24)−(50)(1) c0 = = 1 (2) (2)(−25)−(−52)(1) c1 = =1 (2)

−

14

Para la fila w2 será

(c0 )(50) (c1 )(2) = (c0 ) (c0 )( 52) (0)(2) = (c0 )

d0 = d1 =

Para la fila w1 , e=

− − −

(d0 )(c1 ) (d1 )(c0 ) (d0 )

−

52 52

−

=0

y por fin para la fila w0 , f = d1

quedando la tabla as´ı, 1 w5 4 2 w 3 1 w 52 w2 1 w 0  52 w0

24 50 1 52

− − ≈ −

−25 −52 0

Puede observarse que la pen´ u ltima fila es nula. En este caso se considera que el signo es positivo. Siempre que ocurra que un término de la primera columna es nulo, pero los restantes de su fila sean distintos de cero o inexistentes (como en nuestro ejemplo), se substituye el cero por un n´ umero  > 0, y se siguen calculando los coeficientes de las restantes filas. Cuando se produce este hecho el teorema de Routh-Hurwitz garantiza que el sistema continuo tendr´ a soluciones en el eje imaginario. Por lo tanto el sistema discreto equivalente tendr´ a ra´ıces en la circunferencia unidad (modulo igual a la unidad). Veremos más adelante que esto significa que el sistema tiene un comportamiento oscilatorio no amortiguado. En este ejemplo se ven tres cambios de signo en la primera columna: entre 2 y 1, entre 1 y 52, y entre  y 52. Por lo tanto el sistema tendr´ a tres ra´ıces con parte real positiva.

−

−

−

La aplicación del teorema de Routh-Hurwitz presenta otros casos particulares que deben conocerse. No obstante no voy a incluirlos aqu´ı por el momento.

15

3.

La transformada

Z

La transformada Z transforma una ecuaci´ on en diferencias en una ecuaci´ on polinó mica de una u ´ nica variable compleja. No voy a hacer un estudio exhaustivo de esta transformaci´ on, tan sólo se˜ nalaré los puntos fundamentales que permitan su utilización en la obtención de la ecuación completa de una ecuación en diferencias lineal de parámetros constantes. Dada una variable o función cualquiera yk , utilizaré la notación Z (yk ) = Y (z ) para representar su transformada Z, donde z es una variable compleja genérica2 : ∞ Y (z ) = Z (yk ) = i=0 yi z −k



La transformada Z de una variable desplazada yk−1 viene dada por Z (yk−1 ) = z −1 Y (z )

mientras que Z (yk+1) = zY (z )

− zy

0

y Z (yk+2 ) = zZ (yk+1 )

− zy

1

= z 2 Y (z )

2

− z y − zy 0

1

Observamos que cuando se calcula la transformada Z ”en adelanto” aparecen las condiciones iniciales, cosa que no ocurre cuando se calculan ”en atraso”. Esto se debe a que se est´ a utilizando la tranformada Z unilateral, y se supone que todas las ecuaciones se han expresado en la forma causal3 . La transformada Z es una transformaci´ on lineal, por lo que la ecuaci´ o n en diferencias yk+2 + 2yk+1 5yk = 3uk

−

se transforma en la ecuaci´ on polinómica (z 2 + 2z

2

− 5)Y (z) − (z y

0

2 Voy a

+ z (2y0 + y1 )) = 3U (z )

utilizar aqu´ı la llamada transformada Z unilateral, puesto que esta es la que resulta u ´ til en la resolución de ecuaciones en diferencias en un contexto temporal en el que se admita el principio de causalidad. Será muy importante expresar las ecuaciones en su forma causal. Si no se admite este principio, será más conveniente utilizar la transformada Z bilateral, de la que no voy a hacer aqu´ı ninguna menci´ on adicional. 3 Puede obtenerse demanera recursiva el polinomio de condiciones inciales de Z (yk+m )

16

El polinomio de condiciones inciales4 será el responsable de la soluci´ o n de la ecuaci´ on homog´ enea , P 0 (z ) = z 2 y0 + z (2y0 + y1 )

Esta solución se obtiene haciendo U (z ) = 0, y resolviendo la ecuación (z 2 + 2 z

2

− 5)Y (z) = z y H

0

+ z (2y0 + y1 )

La relaci´ on (z 2 + 2z

− 5)Y (z) = 3U (z) P

representa al sistema con condiciones iniciales nulas, lo que nos permitirá obtener una soluci´ on en diferencias5 . on particular de la ecuaci´ La solución completa ser´ a la suma de la homogénea y la particular. Es decir Y (z ) = Y H (z ) + Y P (z ). Por comodidad de notaci´ on escribiré Y (z ) quitando el sub´ındice a Y P (z ) ya que en lo que sigue analizaré la soluci´ on particular. Podemos observar también que el polinomio que multiplica a la salida transformada Y (z ) representa el polinomio caracter´ıstico del sistema, P (z ) = z 2 + 2z

−5

a cuyas ra´ıces se les llama polos del sistema, y que como se ha visto en el apartado anterior coinciden con los autovalores de la matriz de estados del sistema A. Podemos tambi´ en observar que si se conoce la entrada, se conocer´ a U (z ), y de aqu´ı que pueda obtenerse la salida Y (z ) obtenida como Y (z ) =

3

− U (z )

z 2 +2z 5

De la misma manera que se puede calcular la transformada Z puede calcularse la transformada Z inversa, que permite obtener una ecuaci´ on en diferencias a partir de una ecuaci´ on polinómica compleja. Utilizaré la notaci´ on 4 Si

la ecuación en diferencias no se hubiese expresado en la forma causal, se hubiesen perdido las condiciones iniciales. 5 Cuando se est´ a buscando simplemente la solución particular, se está imponiendo que las condiciones iniciales sean nulas, por lo que carece de importancia expresar la ecuación en diferencias en su forma causal.

17

en es una transformaci´ on Z −1 (Y (z )) = yk . La transfromada Z inversa tambi´ lineal, propiedad que se utilizar´ a al resolver ecuaciones en diferencias. De acuerdo con estas definiciones, se cumplir´ a que yk = Z −1 ( z 2 +23z −5 U (z ))

Vamos a ver enseguida algunos ejemplos. Supongamos que la entrada uk = u es constante para todo k = 0, 1, 2, . . . , y que ademá s es una funci´ on causal, es decir que uk = 0, k < 0. Puede demostrarse que la transfromada Z de una constante6 es Z (u) =

u

1 z

−

1

−

= zzu −1

Puede observarse un detalle importante, y es que cuando se tiene un polo en z = 1 el comportamiento del sistema ser´ a el de un escalón. Si z = 1 el comportamiento ser´ a una oscilaci´ on permanente (no amortiguada), excepto en z = 1.

||

Es habitual expresar las transformadas en función de z −1 , y no en función de z , pero esto no nos debe llevar a ninguna confusi´ on en cuanto al principio de causalidad. Se trata simplemente de una representaci´ o n de n´ umeros comk plejos. Por ejemplo la transfromada Z de la funci´ on a , con k 0, y causal k (es decir a = 0, para k < 0) suele expresarse en la forma

≥

Z (ak ) =

1 1 az

−

1

−

y que, si resulta conveniente puede utilizarse en la forma Z (ak ) =

z z a

−

sin que cambie el significado sobre la causalidad de la funci´ on ak . Esta función potencial ak , tiene un polo en z = a, por lo que ser´ a estable si 2 3 a a cero si a < 1 y tenderá a a < 1. La secuencia 1, a , a , a . . . tender´ infinito en caso contrario. Esta es otra forma de entender el concepto de estabilidad. La salida de un sistema representar´ a una secuencia temporal

||

{

}

6A

||

las funciones constantes causales, en el sentido dado, se les llama funciones escalón o ”step functions”, en inglés. El escalón unitario suele representarse como 1( k).

18

obtenida variando el par´ ametro k. Si esta secuencia converge al punto de equilibrio, podrá decirse que la salida es estable. Si en nuestro ejemplo la entrada fuese constante, con u = 1, la salida podrá obtenerse calculando la transformada Z inversa del cociente de dos polinomios en z : 3z Y (z ) = (z 2 +2z − 5)(z −1) A aquellos valores de z que anulan el numerador de Y (z ) se les llama ceros. En el ejemplo anterior hay un cero en el origen. Para el c´ a clculo de los polos y los ceros debe expresarse la funció n en su forma causal. Veamos unos ejemplos. En todos ellos supondremos condiciones iniciales nulas y una función escalón unitaria de entrada. 1. Ning´ un cero y cuatro polos yk+2 + 2 yk+1 Y (z ) =

− 5y

k

= 3uk−2

3

z (z 2 +2z 5)(z 1)

−

−

2. Dos ceros y cuatro polos yk+2 + 2yk+1

− 5y

k

Y (z ) =

= 2uk + 3uk−2

2z 2 +3 5)(z 1)

z (z 2 +2z

−

−

Para que el sistema sea causal el n´ umero de ceros siempre debe ser menor o igual que el n´ umero de polos. Por u ´ ltimo, es u ´ til conocer el teorema del valor final, ya que permite saber cuál es el comportamiento en régimen permanente sin necesidad de resolver ´ la ecuación en diferencias, tan sólo conociendo Y (z ). Este afirma que si yk es estable se cumple que y∞ = limk→∞ yk = limz→1 (1

19

1

− z− )Y (z)

4.

Resoluci´ on de la ecuaci´ on en diferencias utilizando el m´ etodo de la descomposici´ on o expansi´ on en fracciones parciales.

Antes de calcular la transformada Z inversa deberán eliminarse los polos en el or´ıgen. Para ello se puede hacer lo siguiente. Por ejemplo, si Y (z ) =

2z 2 +3 5)(z 1)

z (z 2 +2z

−

−

entonces se puede hacer 2z 2 +3 5)(z 1)

X (z ) = zY (z ) =

(z 2 +2z

−

−

ya que una vez que se haya calculado x(k), se podr´ a deshacer la transformación, teniendo en cuenta que yk = xk−1 ,

k = 1, 2 . . .

yk = 0,

k

≤0

Realizar este paso es imprescindible para resolver el problema de inversi´ on utilizando el método de expansión en fracciones parciales. El desarrollo en fracciones parciales, cuando los polos son simples , consiste en expresar la funci´ on Y (zz ) , en la forma general Y (z ) z

=

n ci i=1 z zi



−

donde n es el orden de la ecuaci´ on en diferencias, y ci son coeficientes constantes que se calculan utilizando la fórmula: ci = [(z

−z) i

Y (z ) z =zi z

]|

donde zi son los diferentes polos de la expansi´ on. Puesto que no puede haber polos en el or´ıgen en berá tener al menos un cero en el or´ıgen.

20

Y (z ) z

, entonces Y (z ) de-

Una vez calculados los coeficientes podr´ a aplicarse la transformada Z inversa directamente a cada uno de los términos de la expansi´ on, ya que Y (z ) =

luego, yk =

n i=1



n ci i=1 1 zi z 1



−

ci zik ,

yk = 0,

−

k = 0, 1 . . . k<0

El ejemplo anterior no puede desarrollarse en fracciones simples ya que no tiene ning´ un cero en el or´ıgen. No obstante puede aplicarse el principio de superposición, haciendo el truco siguiente: Y (z ) =

Lamando

2z (z 2 +2z 5)(z 1)

−

Y 1 (z ) = Y 2 (z ) =

3z − + z2(z2+2z−5)(z−1)

2z (z 2 +2z 5)(z 1)

−

−

3z z 2 (z 2 +2z 5)(z 1)

−

−

y haciendo el cambio X (z ) = z 2 Y 2 (z ) =

3z (z 2 +2z 5)(z 1)

−

−

puede resolverse la ecuaci´ on aplicando el principio de superposici´ on, obteniendo yk = y1k + y2k = y1k + xk−2 , k = 2, 3 . . . yk = 0,

k

≤1

Para el calculo de y1k y xk puede aplicarse el método de descomposici´ on en Y 1 (z ) X (z ) fracciones parciales de z y de z . Podemos observar que X (z ) = 32 Y 1 (z ), por lo que será suficiente con calcular y1k ya que xk = 32 y1k . Por lo tanto yk = y1k + 32 y1(k−2) , yk = 0,

21

k

k = 2, 3 . . .

≤1

Este ejemplo también puede resolverse sin utilizar el principio de superposición, haciendo z (2z 2 +3) Y (z ) = z 2 (z 2 +2z −5)(z−1) y el cambio X (z ) = z 2Y (z ), ya que X (z ) z

=

2z 2 +3 5)(z 1)

(z 2 +2z

−

−

puede desarrollarse en fracciones simples calculando los coeficientes como se ha hecho antes. Este método también es aplicable al caso de polos m´ ultiples, pero la expansión no puede hacerse igual que en el caso de polos simples. Veremos un ejemplo más adelante.

4.1.

Caso de polos reales simples, con al menos un cero en el or´ ıgen

Como ya se ha mencionado, para el cálculo de la solución particular de una ecuación en diferencias podemos imponer que las condiciones iniciales sean nulas. Sea el sistema 8yk + 2yk−1

= 16uk−2

−y −

k 2

La transformada Z (con condiciones iniciales nulas) será, (8 + 2z −1

2

2

− z− )Y (z) = 16z− U (z)

que si se desea, puede escribirse en la siguiente forma, sin más que multiplicar ambos miembros de la igualdad por z 2 , (8z 2 + 2z

− 1)Y (z) = 16U (z)

Supongamos que la entrada es un escal´ on unidad, es decir que 1 1 z

U (z ) =

−

1

−

Entonces Y (z ) =

(8+2z

1

−

22

16z

−z

2

−

2 )(1

−

−z

1)

−

Hay diversas formas de calcular la transformada Z inversa. En este ejemplo utilizaré el método de expansi´ on en fracciones simples. Para ello debe recordarse que todo polinomio puede descomponerse en el producto de factores de sus ra´ıces. Por ejemplo el polinomio P (z ) = 8z 2 + 2 z

Por lo tanto Y (z ) =

1 )(z 2

− 1 = 8(z +

(1+ 12 z

2z

1 ) 4

2

−

− 14 z

1 )(1

−

−

1 )(1

−z

−

1)

−

Esta expresi´ on puede a su vez escribirse en la forma de una expansión en fracciones parciales, Y (z ) =

c1

1+ 12 z

1

−

+ 1− c12z 4

1

−

+ 1−cz3

1

−

donde c1 , c2 y c3 son coeficientes constantes que pueden obtenerse utilizando la relación ci = [(1 zi z −1 )Y (z )]|z =zi

−

donde zi son los diferentes polos de la expansi´ on. Por lo tanto, para z1 = −21 c1 = [(1

1

2

= [ (1− 1 z 2z1 )(1−z −

− z z− )Y (z)]| 1

z =z1

es decir c1 =

(1

−

4

2z1−2 1 −1 z )(1 4 1

−z1 1 ) = −

−

1)

−

]|z =z1

16 9

Para z2 = 14 , se obtiene c2 =

y para z3 = 1,

2z2 2 (1+ 12 z2 1 )(1 z2 1 ) −

−

c3 =

−

= −932

−

2 (1+ 12 )(1

− 14 ) =

16 9

Una vez descompuesta la ecuación en fracciones parciales la transformada Z inversa es trivial. Esto es as´ı porque, como se ha se˜ nalado en el apartado k anterior la transfromada Z de a es Z (ak ) =

23

1 1 az

−

1

−

y porque puede hacerse uso de la propiedad de linealidad de la transformada Z inversa. De esta manera se obtiene la soluci´ on particular de la ecuación en diferencias, yk = c1 ( −21 )k + c2 ( 14 )k + c3 , yk = 0,

k = 0, 1, 2, . . .

k<0

Podemos analizar el comportamiento de r´ egimen transitorio de la salida observando la evoluci´ on en el tiempo del efecto de cada uno de los polos del sistema. Vemos que el polo z = −21 da una componente de la salida que es oscilatoria amortiguada, mientras que el polo z = 14 da una componente exponencial amortiguada. En ambos casos los polos son estables ya que su módulo es menor que la unidad, pero el comportamiento de cada componente es cualitativamente distinto. Por u´ltimo, el polo en z = 1 produce un escalón constante. Esta superposici´ on de comportamientos nos permite asegurar que cuando k , la salida tenderá a c3 .

→∞

Debe quedar claro que la soluci´ on obtenida es una solución particular y no ´ la solución completa. Esta se obtendr´ a sumando la soluci´ o n general de la ecuación homogénea. Puede comprobarse que se han supuesto condiciones iniciales nulas haciendo k = 0 en la solución: y0 = c1 + c2 + c3 = 0. También puede comprobarse que para k = 1, y(1) = c21 + c42 + c3 = 0.

−

Otro estudio que conviene hacer es el comportamiento en r´ egimen per. Ya se ha hecho este estudio después de manente, es decir cuando k calcular la salida yk , pero puede hacerse antes de invertir Y (z ) aplicando el teorema del valor final. Aplicando este teorema se ve que

→∞

y∞ = limz→1 (1

1

− z− )

2z

(1+ 12 z −1 )(1

−

2 1 −1 z )(1 4

−

−z

1)

−

=

16 9

que como era de esperar coincide con c3 . Este ejemlplo ha servido para exponer un método general basado en la transformada Z de análisis de la salida o solución particular de una ecuación en diferencias, cuando los polos son reales y distintos. Se ha visto as´ı mismo que si éstos son negativos se producir´ a un comportamiento transitorio oscilatorio amortiguado, y si son positivos un comportamiento amortiguado. Tambi´ en se ha visto el comportamiento en régimen permanente. 24

La salida transitoria es una superposici´ on de funciones dependientes de cada uno de los polos, y aunque en este ejemplo no se aprecia, es interesante observar que el efecto de cada una de ellas depende de la rapidez con la que se amortigua la se˜ nal. Puede haber componentes m´ as dominantes que otras. Por ejemplo si un polo es muy peque˜ no comparado con otro, su contribució n a la salida ser´ a peque˜ na ya que se amortiguar´ a rápidamente. Cuanto mayores sean los polos (m´ as proximos en módulo a la unidad) más tarde se amortiguará su contribución, y en consecuencia, m´ as tarde se alcanzar´ a el régimen permanente. Para terminar, no debe olvidarse que este ejemplo se ha resuelto suponiendo que la entrada es un escal´ on unitario. Si la entrada fuese distinta el resultado será, en consecuencia, distinto, pero el efecto de los polos del sistema será el mismo como puede apreciarse de la expansi´ on en fracciones parciales. Las entradas a˜ naden nuevos polos que, como se ha visto en este ejemplo, afectan significativamente al comportamiento en régimen permanente. Hemos visto que si los polos son negativos se producir´ a una oscilaci´ o n (de hecho, una respuesta alternante). Esto puede expresarse a través de la transformada Z inversa como 1 ak cos(kπ ) = Z −1 [ 1+az

1

−

]

donde se supone que a > 0.

4.2.

Soluci´ on con polos complejos simples

Cuando los polos son complejos, éstos aparecen en pares conjugados. Por lo tanto solo puede haber polos complejos si los sistemas son de orden superior o igual al segundo. El polinomio caracter´ıstico de un sistema de segundo orden con ra´ıces complejas conjugadas es: P (z ) = z 2

− 2σz + (σ

2

+ ω2)

siendo las ra´ıces de la ecuaci´ on caracter´ıstica z1 = σ + iω y z2 = z1∗ = σ El sistema será estable si

√σ

2

+ ω2 < 1.

Puede comprobarse que P (z ) = (z

2

− z )(z − z∗) = (z − σ) 1

1

25

+ ω2

− iω.

Para entender los resultados que expondré a continuaci´ on conviene recordar la fórmula de Euler, que representa la ecuaci´ on de una circunferencia de radio unidad: eiφ = cos(φ) + i sin(φ), φ R

∀ ∈

y también que si x = σ + i ω

∈ C, podemos escribir x = |x|e iφ

donde

ω σ

tg (φ) =

|x| = √σ

2

+ ω2

Por u ´ltimo, tambi´ en debe recordarse que xk = x k eikφ

||

Veamos un ejemplo. Sea la ecuación en diferencias yk

−y − − k 1

1 y 2 k 2

− = uk−2

La transformada Z (con condiciones iniciales nulas) será, (1

1

− z−

+ 12 z −2 )Y (z ) = z −2 U (z )

La ecuaci´ on caracter´ıstica es z 2 conjugadas, 12 i 12 .

±

−z+

1 2

= 0, cuyas ra´ıces son complejas

Si suponemos que la entrada es un escalón unidad, U (z ) = 1−1z 1 , por lo que −

Y (z ) =

z 2 −

(1 z

−

1 )(1

−z1z

−

1 )(1

−z1 z

−

∗

1)

−

Para analizar el comportamiento en régimen permanente o estacionario puede aplicarse el teorema del valor final, y∞ = limz →1 (1

1

− z− )Y (z) = lim → z

es decir y∞ =

1 (σ 1)2 +ω 2

−

26

z 2 −

1 (1 z1 z

−

1 )(1

−

=2

−z1 z ∗

1)

−

=

1 (1 z1 )(1 z1 )

−

−

∗

Para el análisis del comportamiento en régimen transitorio, y la obtenci´ on de la solución de la ecuaci´ on en diferencias, puede descomponerse la expresi´ on de Y (z ) en fracciones parciales como c1

Y (z ) =

1 z1 z

−

+ 1−zc2z 1

1

∗

−

1

−

+ 1−cz3

1

−

Ahora los coeficientes ci será n n´ umeros complejos, debido a que z1 es complejo. Vamos a comprobar a continuaci´ on que los dos primeros coeficientes son n´ umeros complejos conjugados c2 = c∗1, y que el tercero, c3 es un n´ umero real. Podemos calcular estos coeficientes aplicando la f´ ormula que ya hemos visto para el caso de ra´ıces reales simples, ci = [(1

1

− z z− )Y (z)]| i

z =zi

donde zi son los diferentes polos de la expansi´ on. c3 = [ (1−z1 z

z 2 −

1 )(1

−

−z1 z ∗

c3 =

1)

−

] z=1 =

|

1 (1 z1 )(1 z1 )

−

−

∗

1 (σ 1)2 +ω 2

−

que es un n´ umero real. En nuestro ejemplo, c3 = 2. 1

z =z1

=

1 (z1 1)(z1 z1 )

z =z1

=

1 (z1 1)(z1 z1 )

− z z− )Y (z)]| = [(1 − z ∗ z − )Y (z )]|

c1 = [(1

1

c2

1

1

de lo que resulta que c2 = c∗1 =

∗

−

∗

−

−

∗

∗

−

1 (i (σ 1)+ω )2ω

− −

Multiplicando numerador y denominador por el conjugado del denominador, se obtiene i (σ −1)−ω c2 = c∗1 = 2ω(ω2 +(σ−1)2 ) que pueden ser expresados en su parte real e imaginaria, cR 1 = cI 1 =

En nuestro ejemplo, cR 1 =

1 2(ω 2 +(σ 1)2 )

−

−

1 σ 2ω (ω 2 +(σ 1)2 )

−

−

−c = −1, y |c | = √2. I

1

1

27

La transfromada Z inversa de Y (z ) será, yk = c1 z1k + c2 (z1∗ )k + c3

que deber´ a ser necesariamente un n´ umero real. Puesto que z1 = z1∗ yk

| | | | = c |z | e + c∗ |z | e− 1

1

k ikφ

1

1

k

ikφ

+ c3 = z1 k (c1 eikφ + c∗1e−ikφ ) + c3

| |

Puede observarse que salvo el término constante c3 , la solución de la ecuaci´ on consiste en una oscilación amortiguada, ya que c1 eikφ + c∗1 e−ikφ = ((c1 + c∗1 )cos(kφ ) + i(c1

− c∗)sin(kφ) = 2c

R

1

1

cos(kφ )

I

− 2c sin(kφ) 1

I donde c1 = cR 1 + i c1 .

Escribiendo c1 en la forma c1 = c1 ei δ

| |

podemos escribir la anterior expresi´ on como c1 eikφ + c∗1 e−ikφ = 2 c1 cos(kφ + δ )

En nuestro ejemplo, δ =

π

2

+

π

4

| | √ , |c | = 2. 1

Por lo tanto la solución de la ecuaci´ on en diferencias será, yk = 2 c1 z1 k cos(kφ + δ ) + c3 ,

| || |

yk = 0,

k = 0, 1 . . .

k<0

Puesto que la funci´ on coseno siempre est´ a acotada, la estabilidad dependerá del valor de z1 . Si éste es menor que la unidad el sistema ser´ a estable, e y∞ = c3 . Si z1 > 1 el sistema ser´ a inestable con oscilaciones de amplitud creciente alrededor de c3 . Y si, por fin, z1 = 1 entonces se producir´ a una oscilación permanente de amplitud constante alrededor de c3. Hay dos casos paticulares que ser´ an analizados más adelante: cuando z1 = 1, en cuyo caso habrá tres polos iguales, y cuando z1 = 1, en cuyo caso habr´ a dos polos iguales.

| |

| |

| | −

28

En nuestro ejemplo yk =

−2√2|z | sin( 1

pi

k

4

yk = 0,

con z1 =

| |

√2 2

pi

k+

4

) + 2,

k = 0, 1 . . .

k<0

< 1.

Puede comprobarse que para k = 0, se cumple que y (0) = 0, condición que hab´ıamos impuesto. También para k = 1 puede verse que y (1) = 2 2 z1 + 2 = 0.

−√| |

Puede demostrarse que se cumplen las dos relaciones generales siguientes, φ R: sin(φ)z 1 Z (sin(kφ )) = 1−2cos(φ)z 1 +z 2

∀ ∈

−

−

Z (cos(kφ )) =

4.3.

−

1 cos(φ)z 1 1 2cos(φ)z 1 +z

−

−

−

−

2

−

Soluci´ on con polos dobles

Hasta ahora sólo se ha hecho un estudio cuando los polos son simples, tanto reales como complejos. En este caso el método de la transformada Z ha sido muy u ´ til. La salida Y (z ) se representa f´ acilmente como una descomposici´ on o expansión en fracciones parciales (siempre que puedan calcularse los polos) cuyos coeficientes se pueden obtener mediante la relaci´ on ci = [(1

1

− z z− )Y (z)]| i

z =zi

donde zi son los diferentes polos de la expansi´ on. La aplicación de la transformada Z inversa permite obtener, de una manera muy sencilla, la solución de la ecuaci´ on en diferencias yk . Sin embargo la relación anterior deja de ser v´ a lida cuando los polos son múltiples, tanto si son reales como complejos. Veamos un ejemplo de polos dobles reales. El caso de polos complejos se resuelve de la misma forma, con la particularidad de que siempre aparecen en pares conjugados. 29

Sea la ecuación en diferencias siguiente:

−y −

yk

k 1

+ 14 yk−2 = uk−2

La transformada Z (con condiciones iniciales nulas) será, (1

1

− z−

+ 14 z −2 )Y (z ) = z −2 U (z )

La ecuaci´ on caracter´ıstica es z 2 z + 14 = 0, cuyas ra´ıces son reales con un orden de multiplicidad igual a 2: z1 = 12 doble.

−

La descomposici´ on en fracciones parciales (cuando z1 = 1 y U (z ) es una función escalón unidad) adopta la forma:



Y (z ) =

c1 1 z1 z 1

−

−

+ (1−cz21zz

1

−

1 )2

−

+ 1−cz3

1

−

Para aplicar la transformada Z inversa es necesario saber que la transformada on temporal ”rampa” causal7 es: Z de una funci´ Z (k) =

z 1 (1 z 1 )2 −

−

−

y que

z 1

Z (ka k−1 ) =

−

(1 az

−

1 )2

−

Por lo tanto, aplicando la transformada Z inversa se obtiene la solució n de la ecuación en diferencias (con condiciones iniciales nulas): yk = c1 z1k + c2 kz 1k−1 + c3 , yk = 0,

k = 0, 1 . . .

k<0

Podemos observar que aparece un término adicional con respecto al caso de polos simples, de la forma kz1k−1 , lo cual parece que puede afectar a la estabilidad del sistema. Se puede analizar esta cuesti´ on de la siguiente manera. 7 La

rampa causal se define como la función temporal yk = k, k 0, e yk = 0, para abola” causal y k < 0. As´ımismo pueden definirse otras funciones causales como la ”par´ siguientes potencias de k. La transformada Z de una función potencia m-sima de k es

≥

Z (km ) =

1

Q(z (1−z 1 ) z

−

1

−

−

)

+1

m

donde Q(z −1 ) es un polinomio en z −1 , que para m = 2 es Q(z −1 ) = 1+ z −1 , y para m = 3 es Q(z −1 ) = 1 + 4z −1 + z −2 .

30

Como kz1k−1 = z11 (kz1k ), cuando z1 < 1 el sistema será estable, ya que el crecimiento de la función k es más lento que el decrecimiento de z1k , es decir que limk→∞ (kz1k ) = 0 cuando z1 < 1. Por lo tanto los polos m´ ultiples no afectan a las condiciones de estabilidad, pero obviamente s´ı afectan al comportamiento en régimen transitorio: la curva ya no es una exponencial decreciente.

| | | |

El cálculo de los coeficientes para el caso de polos dobles se hace de la siguiente forma: Y (z ) c2 = [(1 z1 z −1 )2 z 1 ]|z=z1 d c1 = [ dz

En el ejemplo

− {(1 − z z− )

d c1 = [ dz

1 2 Y (z )

1

c2 =

z1 1 −

=

1 z1 1 −

−

z 1 1 z 1 −

{ − }]| −

z =z1

z 1 −

}]|

z =z1

1 z1 1

−

1

d = [ dz

{ − }]| z 1

z =z1

1 (z1 1)2

− −

c1 =

que en nuestro ejemplo es c2 =

−

−2 y c = −4. 1

Para el cálculo de c3 se utiliza la fórmula de polos simples, y en nuestro ejemplo resulta c3 = (1−11 )2 = 4. 2

Por lo tanto la solución de la ecuaci´ on en diferencias del ejemplo es yk =

1 k 2

−4( ) − 2k(

1 k 1 ) 2

− + 4, k = 0, 1 . . .

yk = 0,

k<0

que puede escribirse en la forma: yk =

−4(1 + k)2−

k

+ 4,

yk = 0,

k = 0, 1 . . .

k<0

Podemos comprobar que y (0) = c1 + c3 = 0, como hab´ıamos impuesto. También para k = 1, se cumple que y(1) = 0. La descomposici´ on en fracciones parciales y el cálculo de sus coeficientes para el caso de polos con un orden de multiplidad mayor que dos es engorroso, por 31

lo que conviene utilizar alguna otra técnica de inversi´ on de Y (z ). La forma más general de resolver la ecuaci´ on en diferencias utilizando la transfrormada Z es utilizar directamente la transformada Z inversa a través de la integral de inversión. No voy a entrar en detalles sobre esta cuesti´ on, tan sólo explicaré el método en el apartado siguiente.

5.

Resoluci´ on de la ecuaci´ on en diferencias utilizando el teorema de los residuos

Una vez que se ha calculado Y (z ) puede obtenerse yk como yk =



n i=1

ri

donde n es el n´ umero de polos distintos de la ecuaci´ on en diferencias (si todos los polos son simples, coincidir´ a con el orden del sistema) y ri son los residuos k−1 de Y (z )z en cada uno de los polos zi , i = 1 . . . n. Los residuos se calculan de la siguiente forma: ri =

1

limz →zi (m−1)!

{

dm 1 dz m 1 −

−

[(z

m

−z) i

Y (z )z k−1 ]

donde m es el orden de multiplicidad del polo zi , y (m (m 1) (recuérdese que 1! = 0! = 1).

−

}

− 1)! es el factorial de

1. Polos dobles En el ejemplo anterior, la ecuación en diferencias era: yk

+ 14 yk−2 = uk−2

−y −

k 1

cuya transformada Z (con condiciones iniciales nulas, y una entrada escalón unitario) es, Y (z ) =

z

(z 2

−z+ 14 )(z−1)

Vemos que hay un polo doble en z1 = Por lo tanto, n = 2 y k 1

1 2

y un polo simple en z2 = 1.

{ − 1)Y (z)z − } = lim → {

r2 = limz→1 (z

32

z

z

1 z 2 z+ 1 z 4

−

k 1

−

}=4

r1 = limz→ 1 2

r1 = limz → 1 2

{

r1 = limz→ 1 2

{

d dz

d zk dz z 1

{

[(z

−

1 2 ) Y (z )z k 1 ] 2

−

[ − ] = limz→ 1 2

}

z k (kz 1 (z 1) 1) (z 1)2 −

−

− −

}

kz k 1 (z 1) z k (z 1)2 −

{

− − −

}

k

} = −4(1 + k)2−

luego yk = r1 + r2 , que coincide con lo obtenido en el apartado anterior. 2. Polos complejos Sea la ecuaci´ on en diferencias + 12 yk−2 = uk−2

−y −

yk

k 1

La transformada Z (con condiciones iniciales nulas y entrada escal´ on unidad) es: Y (z ) = (z 2 −z +z1 )(z −1) 2

que tiene dos ra´ıces complejas conjugadas simples, z1 = 12 + i 12 y z2 = z1∗ = 12 i 12 , y una ra´ız real simple z3 = 1, por lo que n = 3, con los residuos

−

k 1

r1

z

k 1

{ − 1)Y (z)z − } = lim → { − z − } = 2 = lim → {(z − ( + i ))Y (z)z − }

r3 = limz→1 (z

z

( 12 +i

r1 = limz →z1

1 z 2 z+ 1 2

z

1 2

1 ) 2

1 2

k 1

zk (z 1)(z z1 )

z1k (z1 1)(z1 z1 )

zk (z 1)(z z1 )

(z1

− − }= − − r = lim → − − }= − − Haciendo z = |z |e = σ + i ω y z − 1 = ae , donde a = |z − 1|, y teniendo en cuenta que |z | = |z ∗ | y |z − 1| = |z ∗ − 1|, podemos 2

1

1

z

z1

∗

{ {

∗

iφ

∗

∗

iα

1

1

(z1 )k 1)(z1 z1 ) ∗

∗

1

1

1

1

escribir los residuos relacionados con los polos complejos conjugados en la forma |z1|k eikφ = |z1 |k ei(kφ α) r1 = 2iωae iα aω 2i −

r2 =

por lo que r2 = r1∗ , y

z1 k e 2iωae

−|

|

r1 + r2 =

ikφ

−

iα

−

=

−|

z1 k e aω

|

i(kφ−α)

−

2i

|z1 |k sin(kφ − α) |z1−1|ω 33

Puede comprobarse que la soluci´ on yk = r1 + r2 + r3 coincide con la obtenida anteriormente en el apartado 4.2.8

8 Como

c1 =

1 (z1 −1)(z1 −z1 ) ∗

= c1 eiδ , entonces c1 =

| |

| |

34

1 , 2|z −1|ω

yδ=

π

− 2 − α.

Ecuaciones en Diferencias

Recommend Documents