Control Digital

CONTROL DIGITAL DOCENTE: AGUSTÍN SOTO OTALORA INGENIERO ELECTRÓNICO ESP. AUTOMATIZACIÓN INDUSTRIAL MAGISTER EN INGENIER

Views 200 Downloads 2 File size 2MB

Report DMCA / Copyright

DOWNLOAD FILE

Recommend stories

Citation preview

CONTROL DIGITAL

DOCENTE: AGUSTÍN SOTO OTALORA INGENIERO ELECTRÓNICO ESP. AUTOMATIZACIÓN INDUSTRIAL MAGISTER EN INGENIERÍA DE CONTROL INDUSTRIAL

UNIVERSIDAD SURCOLOMBIANA FACULTAD DE INGENIERÍA PROGRAMA INGENIERÍA ELECTRÓNICA NEIVA - 2013

PRÓLOGO Este módulo de control digital no es un libro, es la recopilación de apuntes que durante 5 años que he venido dictando esta asignatura como docente en el programa de Ingeniería Electrónica de la Universidad Surcolombiana. Para ello se han tomado como guía muchos libros que aparecen como referencia al final de cada capítulo. Este módulo pretende suplir las necesidades de nuestros estudiantes, que en su mayoría son de estratos 1 y 2 y no tienen la capacidad económica para adquirir un ejemplar, aunque la biblioteca de la Universidad tiene el número que recomienda el ministerio de educación nacional, no es suficiente. El objetivo que se ha perseguido al redactar este módulo es que el estudiante tenga a mano desarrollado todo el contenido curricular del curso CONTROL DIGITAL, cosa que es difícil encontrar en un solo libro, aunque en la actualidad hay muy buenos textos. El módulo de control digital intenta exponer, en forma razonada y clara la teoría y práctica del control digital al emplear en el control de procesos un computador o un microprocesador en los lazos de control. Los ochos capítulos desarrollados en este módulo tienen la particularidad que no tiene ningún libro hasta ahora y es que se hace un desarrollo tanto analítico y por mat-lab de todos los problemas de análisis como de diseño.

AGUSTÍN SOTO OTÁLORA

CAPÍTULO I

3

AGUSTÍN SOTO OTÁLORA

1. ECUACIONES EN DIFERENCIA Y TRANSFORMADA Z 1.1 ECUACIONES EN DIFERENCIA Una ecuación en diferencias es una ecuación que relaciona varias secuencias con ellas mismas desplazadas (o sea, una ecuación de recurrencia).

a 0 y k  a1 y k 1  a 2 y k  2    a n y k  n  b0 u k  b1u k 1  b2 u k  2    bm u k  m n

a y i 0

i

m

k i

  b j uk j j 0

1.1.1 DIFERENCIA FINITA

Sea

y  f [k ]

una función para k  Z, se llama primera diferencia, o diferencia finita de

primer grado de y  f [k ] a la expresión dada por

 f [k ]  f [k  1]  f [k ] , gráficamente es

FIGURA 1.1

De esta manera vemos que

 f [ k ] corresponde al incremento que sufre y  f [k ] cuando la

variable k se incrementa en una unidad.

CONTROL DIGITAL

6

Ejemplo: Si

y  [k ]  1  2k  3k 2 hallar la diferencia de primer orden. Solución:

Yk  f [k ]  1  2k  3k 2 Yk 1  f [k  1]  1  2(k  1)  3(k  1) 2  1  2k  2  3(k 2  2k  1)  6  8k  3k 2

 f [k ]  f [k  1]  f [k ]  (6  8k  3k 2 )  (1  2k  3k 2 )  5  6k

Como se observa, a partir del resultado del ejemplo anterior,  f [ k ] es otra función que depende de la misma variable k, por tanto se puede hablar de la segunda diferencia de f [k ] (o primera diferencia de  f [ k ] ). Esta se escribe:

 2 f [k ]   f [k  1]   f [k ] y así sucesivamente para diferencias de orden mayor.

EJERCICIOS Hallar la primera y segunda diferencia de cada una de las siguientes funciones.

1.

f [k ]  k 2  2k  1

2.

g [k ]  1 / k

3.

h [k ]  ln (k )

UNIVERSIDAD SURCOLOMBIANA

CAPÍTULO I

7

1.1.2. ECUACIONES DE DIFERENCIA FINITA DE PRIMER ORDEN Una ecuación que relacione los valores de una función

y  f [k ] con una o varias de sus

diferencias finitas, se llama ecuación de diferencia finita en f [k ] . En adelante se llamará simplemente ecuaciones de diferencia. El orden de estas ecuaciones lo determina la diferencia de mayor grado que se encuentre en la ecuación.

Ejemplos: 1.

3Yk 1  Yk  0

2.

Yk 1  4Yk  k

3.

Yk 1  Yk  3 (2) k

1.1.3. MODELAJE DE SISTEMAS DE TIEMPO DISCRETO MEDIANTE ECUACIONES EN DIFERENCIAS Como los sistemas de tiempo discreto que se usan aquí no son tomados del mundo físico real (la naturaleza) sino que son una invención del hombre, normalmente no parten de un modelo físico, como ocurría en los sistemas de tiempo continuo, sino que directamente se representan mediante algún tipo de modelo matemático.

Lo más normal es que los sistemas de tiempo discreto se representen gráficamente mediante diagramas de bloques, donde los elementos fundamentales son: multiplicación por una constante, sumador y retardo. Así, un sistema de tiempo discreto puede representarse mediante un diagrama de bloques con estos elementos básicos y de él se puede obtener la ecuación en diferencias correspondiente. Se asume que uk es la entrada, yk es la salida y n es el orden del sistema de tiempo discreto.

AGUSTÍN SOTO OTÁLORA

CONTROL DIGITAL

8

Ejemplo: Obtener un modelo matemático para los voltajes de nodo Vk (k= 1, 2, 3, …, N) en la red de resistencias R-2R.

FIGURA 1.2 Solución: Planteando la ecuación de suma de corrientes en l nodo k-1 (voltajes de nodo) se obtiene:

1 1 1 1 1  R  2 R  R V k 1  R V k  2  R V k  0  

5 1 1 V k 1  V k  2  V k  0 2R R R



O mejor

2Vk  5Vk 1  2Vk 2  0

para k=2, 3, 4, …, N

Que es la ecuación en diferencias que modela a los voltajes de nodo para la red R-2R. En el nodo N se tiene la restricción impuesta por la ecuación

VN 

1 V N 1 2

Ejemplo: Obtener la ecuación en diferencias que modela al sistema de tiempo discreto de la siguiente figura:

FIGURA 1.3

UNIVERSIDAD SURCOLOMBIANA

CAPÍTULO I

9

Solución: Del diagrama de bloques se obtiene:

y k  u k  ay k 1

Y reorganizando queda la ecuación de diferencias así:

y k  ay k 1  u k

Ejemplo: Obtener la ecuación en diferencias que permite modelar al sistema de tiempo discreto representado por el diagrama de bloques de la siguiente figura:1.4

FIGURA 1.4

Solución:

Del diagrama de bloques se tiene que:

yk 

quedando la ecuación de diferencias así:

yk 

AGUSTÍN SOTO OTÁLORA

1 2 y k 2  y k 3  u k  2u k  2 4 8

1 2 y k 2  y k 3  u k  2u k  2 4 8

CONTROL DIGITAL

10

1.1.4. SOLUCIÓN DE ECUACIONES DE DIFERENCIA FINITA Una función Y  f [k ] es una solución de la ecuación de diferencia finita si está definida para todo k= 0, 1, 2, 3,... y satisface la ecuación dada.

Ejemplo: Probar que

Yk  f [k ]  k 2  1

,

.

es una solución de la ecuación de diferencia:

Yk 1  Yk  2k  1 Solución:

Yk 1  Yk  [( k  1) 2  1]  [k 2  1]  2k  1 o sea que la función Y  f [k ] cumple con la ecuación de diferencia dada en el ejemplo.

Existen dos clases de soluciones para una ecuación de diferencia finita: la solución particular y la solución general. Se llama solución particular a aquella que no tiene constantes arbitrarias y cumple con la definición anterior de solución, como en el caso del ejemplo. Se llama solución general a aquella que teniendo una o varias constantes arbitrarias, cumple con la definición anterior de solución. Ejemplo:

Yk  f [k ]  k 2  1  C , donde C es una constante arbitraria, es también una

solución de

Yk 1  Yk  2k  1 .

EJERCICIOS. Dada la función Yk  f [k ] averiguar si es la solución de la ecuación correspondiente 1.

Yk  k

;

Yk 1  Yk  1

UNIVERSIDAD SURCOLOMBIANA

CAPÍTULO I

11

2.

Yk  4 (2) k

;

Yk 1  2Yk  0

3.

Yk  1  (k  1) * k

;

Yk 1  Yk  2k  3

4.

Yk  3

;

Yk 1  Yk  0

5.

Yk 

(1  k ) * k C 2

;

Yk 1  Yk  k

6.

Yk 

1 1 k

;

Yk 1  3Yk  1

7.

8.

9.

1  3k 1 Yk  2  1 Yk  1      2

Y  ; k 1

Yk 1  Yk

k

Yk  3(2k 1  1)

k Y  e k 10.

;

2Yk 1  Yk  3

;

Yk 1  2Yk  1

;

Yk 1  Yk  (1  e)Yk

1.1.5. ECUACIÓN DE DIFERENCIA LINEAL DE PRIMER ORDEN Se llama ecuación de diferencia lineal de primer orden con coeficientes constantes en Yk a toda expresión de la forma:

a1 Yk 1  a0 Yk  f [k ] AGUSTÍN SOTO OTÁLORA

CONTROL DIGITAL

12

con a0, a1 constantes y f [k ] una función tal que k = 0, 1, 2, 3, … Ejemplos: 1.

5Yk 1  3Yk  k  1

2.

2Yk 1  Yk  2k

La ecuación en diferencia de orden n es

an Yk  n  an 1 Yk  n 1    a1 Yk 1  a0 Yk  U k Se recalca que estas ecuaciones son aplicables si y sólo si el valor de la función en un período cualquiera, está relacionado o depende del valor de dicha función en el período inmediatamente anterior, además la variable independiente debe tomar los valores enteros 0, 1, 2, 3, … Hasta aquí siempre se ha considerado el intervalo [k, k+1] y sus valores correspondientes de Yk y Yk+1 en los extremos de este intervalo. Sin embargo, se puede presentar otra notación, considerando el intervalo [k-1, k], cuyo valor de Yk  f [k ] serán Yk-1 y Yk respectivamente, por tanto

 f [k ]  Yk 1  Yk y la ecuación de diferencia lineal de primer orden con coeficientes constantes se escribirá:

a1 Yk  a0 Yk 1  f1 [k ]

Donde a0, a1 son constantes y f1 [k ] es una función que depende de k relacionada con f [k ] . El manejo de estas ecuaciones es idéntico al de las anteriores, simplemente que Yk de la primera ecuación corresponde a Yk-1 de la segunda y Yk+1 de la primera ecuación corresponde a Yk de la segunda. Ejemplo: En un instante cualquiera, una señal es igual a las tres cuartas partes de la recibida en el periodo anterior aumentada en 0.2. Plantear una ecuación de diferencia que relacione dos señales consecutivas.

UNIVERSIDAD SURCOLOMBIANA

CAPÍTULO I

13

Solución:

Yk 1  0.75Yk  0.2

ó también:

Yk  0.75Yk 1  0.2

A manera de ejemplo, consideremos una interesante población de caníbales a la cual se le toma muestras de la población cada 15 años, esta crece la mitad en este periodo. Además el muestreo ha demostrado que dichos caníbales siempre se comen (literalmente) el 80% de la población que había hace dos periodos de muestreo, es decir, 30 años y que además, caníbales inmigrantes de otros pueblos llegan a esta población. Este es un típico problema de sistemas discretos y queremos modelarlo mediante un ecuación de diferencia, para ello suponemos un periodo de muestreo de 15 años y llamaremos a y[k ] la salida, es decir, el periodo en que queremos conocer la población de caníbales, donde k puede tomar valores de 0, 1, 2, 3,… etc., y u[k ] la asumimos como la entrada, es decir los caníbales inmigrantes de otros pueblos, así la población de caníbales en cualquier periodo k (recordar que un periodo son 15 años) y[k ] , está dada por:

y[k ]  1.5 y[k  1]  0.8 y[k  2]  u[k ] Donde y[ k  1] y y[ k  2] representan la población de caníbales de hace uno y dos periodos respectivamente. Se puede organizar la ecuación de tal forma que:

y[k ]  1.5 y[k  1]  0.8 y[k  2]  u[k ] O también

y[k  2]  1.5 y[k  1]  0.8 y[k ]  u[k  2] Las anteriores ecuaciones de diferencia modelan el sistema discreto propuesto.

1.1.6. SOLUCIÓN DE LAS ECUACIONES DE DIFERENCIA Uso de la tabla Un método para resolver una ecuación de diferencia, no muy práctico pero interesante para conocer la dinámica del sistema es la construcción de una tabla en la cual se van registrando las salidas y entradas para los respectivos periodos en el tiempo conforme éste avanza. Para explicar este método, retomemos el ejemplo de la población de caníbales y calculemos cual sería la oblación en 60 años, es decir 4 periodos de muestreo, suponiendo que en un instante AGUSTÍN SOTO OTÁLORA

CONTROL DIGITAL

14

inicial la población existente son los 2 caníbales que cada periodo llegan, inmigrantes de otros pueblos que se toman como entrada al sistema. Conociendo la ecuación de diferencia:

y[k ]  1.5 y[k  1]  0.8 y[k  2]  u[k ] Construimos nuestra tabla así:

k

y[k ]

y[ k  1]

y[ k  2]

u[k ]

0

2

0

0

2

1

5

2

0

2

2

7.9

5

2

2

3

9.85

7.9

5

2

4

10.45

9.85

7.9

2

Tabla 1

De la tabla anterior, se llega a la conclusión de que la población de caníbales en 60 años es de 10.45. Como se puede observar, este método es muy conveniente para mirar la dinámica del sistema, pero que tal si la pregunta fuera cual es la población en 1000 años?. La construcción de la tabla para resolver esta inquietud es muy tediosa, de ahí que mejor se empleen métodos matemáticos para resolver este tipo de ecuaciones de diferencia

UNIVERSIDAD SURCOLOMBIANA

CAPÍTULO I

15

SOLUCIÓN DE ECUACIONES EN DIFERENCIA SOLUCIÓN ECUACIÓN HOMOGÉNEA

FIGURA 1.5 Caso 1: Si todas las raíces ri son diferentes, la solución general es: n

y[k ]   c i ri , k

i 1

k 0

De las condiciones iniciales se determinan las constantes ci: n

y[k ]   c i ri , i 1

k

k= 0, 1, 2, …,n -1

AGUSTÍN SOTO OTÁLORA

CONTROL DIGITAL

 y[0]  1  y[1]   r   1  y[2]      n 1  y[3] r1

16

1 r2  r2

n 1

  c1   c   2      n 1    rn  c n 

  

1 rn 

Si todas las ri son diferentes, el conjunto de ecuaciones resultante de esta matriz siempre se puede resolver. Es posible que algunas raíces así como algunas constantes sean complejas conjugadas, pero para condiciones iniciales reales, la solución siempre será real. Caso 2: en caso de que se presenten raíces reales repetidas con multiplicidad m, se tiene que adicionar a la solución y a la matriz, términos como:

ri , kri , , k m1 ri k

k

k

Ejemplo: Resolver la siguiente ecuación homogénea:

y[k  3]  7 y[k  2]  16 y[k  1]  12 y[k ]  0 Con condiciones iniciales y[0]=1, y[1]=3, y[2]=0. La ecuación característica

r 3  7r 2  16r  12  0 tiene las siguientes raíces:

r1=3, r2 = r3 =2 La solución general está dada por:

y[k ]  c1 3 k  c 2 2 k  kc3 2 k

ó

y[k ]  c1 3 k  (c 2  kc3 )2 k

Las constantes se calculan de:

1 1 1 0  c1  3  3 2 2 c      2  0 9 4 8 c 3 

De donde c1=-8, c2=9 y c3=4.5. Obteniéndose la solución: UNIVERSIDAD SURCOLOMBIANA

CAPÍTULO I

17

y[k ]  (8)3 k  (9  4.5k )2 k

Caso 3: En el caso de tener raíces conjugadas complejas, se puede probar que para condiciones iniciales reales, las constantes ci son también conjugadas complejas. Esto resultará en una solución real.

Por ejemplo si:

c i  R 0 e j 0

ri  Re j ,

ri 1  Re (  j ) , c i 1  R0 e

(  j 0 )

De la solución general: n

y[k ]   ci ri , k

i 1

k 0

Se tiene que:

ci ri  ci 1rik1  R0 Rk e j ( k 0 )  R0 Rk e j ( k 0 ) k

ci ri  ci 1rik1  2R0 R k cos(k  0 ) k

AGUSTÍN SOTO OTÁLORA

CONTROL DIGITAL

18

Ejemplo: Supóngase la ecuación de diferencia de segundo orden:

y[k  2]  2ay[k  1]  y[k ]  0 Esta ecuación tiene como ecuación característica asociada a son:

r1 , r2  a  a 2  1 Lo que implica una solución homogénea de la forma:

y[k ]  c1 r1  c 2 r2 k

k

Así:







k

y[k ]  c1 a  a  1  c 2 a  a  1 2

2



k

Dependiendo de a, la solución puede tener diferentes formas: Si a 1 entonces tenemos raíces reales y se puede aplicar el caso 1.

La ecuación no homogénea La solución de la ecuación general se obtiene mediante la adición de una solución particular a la solución de la ecuación homogénea, dicha solución particular no necesariamente satisface las condiciones iniciales. Una solución particular puede encontrarse con el método de variación de parámetros o el método de parámetros indeterminados. Este último método propone una solución particular semejante al lado derecho de la ecuación pero con un parámetro indeterminado . Algunas de las soluciones propuestas dependiendo del lado derecho de la ecuación se muestran en el siguiente cuadro:

AGUSTÍN SOTO OTÁLORA

CONTROL DIGITAL

20

Lado derecho de la ecuación en un instante k

Solución propuesta

k

kn

ak

k nak

k n a k cos k

1k   0

sin(𝛼𝑘) 𝑜 𝑐𝑜𝑠(𝛼𝑡)

𝐴0 𝑐𝑜𝑠(𝑘𝛼 ) + 𝐴1 𝑠𝑖𝑛(𝑘𝛼 )

Polinomio p(k) de grado m

𝐴0 𝑘 𝑚 + 𝐴1 𝑘 𝑚−1 + ⋯ + 𝐴𝑚

n

 k i 0

 ak

i

i

 n     i k i a k  i 0 

 n     i k i a k cos(k   0 )  i 0 

1[ k ]

 1[ k ]

(-1)K

a(-1)K

Tabla 2

Si a ó

e  j

son raíces de multiplicidad m de la ecuación característica, entonces la

solución propuesta debe ser multiplicada por

km.

Ejemplo: Resolver la siguiente ecuación:

y[k  3]  7 y[k  2]  16 y[k  1]  12 y[k ]  (1) k UNIVERSIDAD SURCOLOMBIANA

CAPÍTULO I

21

Con condiciones iniciales y[0]=1, y[1]=3, y[2]=0. El lado derecho de la ecuación es una serie de potencias y por tanto seleccionamos una solución particular de la misma serie de potencias. Así, se puede seleccionar:

y[k ]   (1) k Sustituyendo:

  (1) k 3  7  (1) k  2  16  (1) k 1  12  (1) k  (1) k Lo que resulta en:

  (1) k (1) 3  7  (1) k (1) 2  16  (1) k (1)1  12  (1) k  (1) k De donde

 36  (1) k  (1) k Por tanto   

1 . 36

Teniendo en cuenta que la solución a la ecuación homogénea ya fue calculada en el ejemplo anterior, la solución general es:

y[k ]  c1 3 k  (c 2  kc3 )2 k 

y[k ] 

1 (1) k 36

ó

1 (1) k  c1 3 k  (c 2  kc3 )2 k 36

Ahora calculamos las constantes ci

 37 / 36  1 1 0  c1  107 / 36  3 2 2 c      2   1 / 36  9 4 8 c3 

AGUSTÍN SOTO OTÁLORA

CONTROL DIGITAL

Resultando

c1  

22

156 279 , , c 2   316 , c 3   36 36 36

Retomando ahora el problema de caníbales y conociendo su ecuación de diferencia dada por:

y[k  2]  1.5 y[k  1]  0.8 y[k ]  u[k  2] Podemos encontrar la solución de la anterior ecuación y de esa forma conocer la población de caníbales en cualquier periodo k sin necesidad de construir una tabla. Para ello, empezamos por la solución a la ecuación homogénea y por tanto necesitamos conocer las raíces de la ecuación característica:

r 2  1.5r  0.8  0 Cuyas raíces son: r1 , r2

 0.75  0.4873 j

R  (0.75) 2  (0.487) 2  0.894

o expresadas de otra forma:

y

rad  0.487    33.14  0.576 s  0.75 

  tan 1  De donde:

r1  (0.894)e j 0.576

r2  (0.894)e  j 0.576

Lo que implica una solución homogénea de la forma:

y[k ]  c1 (0.894) k e ( j 0.576) k  c 2 (0.894) k e (  j 0.576) k Ahora solucionamos la ecuación particular y[k ]   : El lado derecho de la ecuación es una entrada paso de la forma 2(1[k]) y por tanto la solución propuesta tendrá la forma (1[k]), así en la ecuación:

y[k  2]  1.5 y[k  1]  0.8 y[k ]  u[k  2] Se sustituye la solución propuesta y por tanto: De donde

 1.5  0.8  2

  6.666 UNIVERSIDAD SURCOLOMBIANA

CAPÍTULO I

23

Por consiguiente la solución total será:

y[k ]  c1 (0.894) k e ( j 0.576) k  c 2 (0.894) k e (  j 0.576) k  6.666 Teniendo en cuenta las raíces complejas la solución se puede escribir como:

y[k ]  2 R0 (0.894) k cos( k 0.576   0 )  6.666 Ahora con condiciones iniciales y[0]=2, y[1]=5 se tiene que:

2  2R0 cos  0  6.666

5  2R0 (0.894) cos(0.576   0 )  6.666 Así tenemos que:

2R0 cos   4.666

2R0 cos(  0.576)  1.864

De las anteriores ecuaciones se obtiene que:

2R0 cos   4.666

2R0 sen  3.763

Y así:

2R0  5.994

  4.666    3.82  5.994 

  cos 1 

De donde la solución total estará dada por:

y[k ]  5.994(0.894) k cos(k 0.576  3.82)  6.666 1.2 DEFINICIÓN DE LA TRANSF_Z Sea

x(t)

: Señal en tiempo continuo

x(kT) : Señal en tiempo discreto AGUSTÍN SOTO OTÁLORA

CONTROL DIGITAL

24

La transformada z se define de la siguiente manera: 

Z[x(kT)] = X(z) =

 x ( kT ) z

k

k 0

= x(0) + x(T) z -1 + x(2T) z

1.3

-2

+ ··· + x(kT) z

–k

+ ···

TRANSFORMADA Z DE FUNCIONES ELEMENTALES

1.3.1

ESCALÓN UNITARIO kT  0 , y

x(kT) =1,

x(kT) = 0,

kT  0

Reemplazando en la definición 

X ( z )   x(kT ) * z k k 0



X ( z )  1* z k 0

k

;

donde x(kT )  1



 z k  sk  1  z 1  z 2  z 3    z k k 0

Si multiplicamos la ecuación (1) en ambos miembros por Z-1

z 1sk  z 1  z 2  z 3  z 4    z  ( k 1) Ec. (2) Si restamos la ecuación (2) - (1)

z 1sk  sk  1  z  k sk ( z 1  1)  1  z  k

 1  z k sk  1 z 1 UNIVERSIDAD SURCOLOMBIANA

Ec. (1)

CAPÍTULO I

25

11  1  1  1 k  z 1 z 1

lim S k  lim

k 

X ( z) 

z z 1

La anterior relación se obtiene por Matlab de dos formas: Utilizando la definición de la transformada Z >> syms z k

%se crean las variables simbólicas z,k

>> f=symsum(z^(-k),k,0,inf); >> pretty(f)

%se arregla la presentación z ------1 + z

Utilizando el comando ztrans

>> syms z k >> g=1^(k);

%se crean las variables simbólicas z,k %se crea la función a transformar

>> f=ztrans(g); %se calcula la transformada z >> pretty(f)

%se arregla la presentación z ------1 + z

Para graficar la señal escalón unitaria discreta por Matlab, se hace: AGUSTÍN SOTO OTÁLORA

CONTROL DIGITAL

26

% GENERACIÓN DE ESCALÓN UNITARIO DISCRETO x = ones (1,11);

% define once valores de 1's

v = [0 10 0 2];

% define valores de ejes

axis (v); plot (x,'ro')

% grafica círculos de color rojo

xlabel ('k')

% asigna rotulo al eje x

ylabel ('x(k)')

% asigna rotulo al eje y

title (‘ESCALON UNITARIO DISCRETO’)

FIGURA 1.6

1.3.2

RAMPA UNITARIA

% GENERACIÓN DE LA RAMPA UNITARIA DISCRETA k = 0:10;

% define valores de k

x = k;

% función rampa para x UNIVERSIDAD SURCOLOMBIANA

CAPÍTULO I

27

axis([0 10 0 10]);

% define ejes

grid

% rejilla para grafica

plot(k, x,'ro')

% grafica x en función de k

xlabel('k');

% rotulo para eje x

ylabel('x(k)');

% rotulo para eje y

title('RAMPA UNITARIA DISCRETA')

FIGURA 1.7

kT, x (kT)    0,

para k  0, 1, 2, 3, ··· señal dis creta para k  0

X(z)  Z[x(kT)]  Z[kT] 



 kT * z

-k

k 0

AGUSTÍN SOTO OTÁLORA

CONTROL DIGITAL

28

sk  0  z 1  2 z 2  3z 3    kz  k

Ec. (1)

Si multiplicamos la ecuación (1) en ambos miembros por z-1

z 1sk  z 2  2 z 3  3z 4    (k  1) z  ( k 1)

Ec. (2)

Si restamos la ecuación (2) - (1)

z 1sk  sk  z 2  2 z 3  3z 4    (k  1) z  ( k 1)  z 1  2 z 2  3z 3    kz k

sk ( z 1  1)   z 1  z 2  z 3    z  k Del ejercicio anterior se tiene que

sk 

 (1  z 1  z  2  z  3    z  k ) 

z 1 Tz * 1  z  1 z  1 ( z  1)2

X ( z) 

Tz 1

1  z 

1 2



Tz z  12

La anterior relación se obtiene por Matlab de dos formas: Utilizando la definición de la transformada Z >> syms z k

%se crean las variables simbólicas z,k

>> f=symsum(z^(-k)*k,k,0,inf); >> pretty(f)

%se arregla la presentación z --------(-1 + z)2

Utilizando el comando ztrans >> syms z k

%se crean las variables simbólicas z,k UNIVERSIDAD SURCOLOMBIANA

z z 1

CAPÍTULO I

29

>> g=k;

%se crea la función a transformar

>> f=ztrans(g); %se calcula la transformada z >> pretty(f)

%se arregla la presentación z --------(-1+z)2

1.3.3

POTENCIAL: a

k

(a = constante)

%GENERACION DE LA FUNCION POTENCIAL

x(k) = 2

k

k=linspace(0,5,20);

% define valores de k

x=2.^ k;

% función potencial

grid

% rejilla para gráfica

plot(k, x,'ro')

% gráfica x en función de k

xlabel('k');

% rotulo para eje x

ylabel('x(k)');

% rotulo para eje y

title('POTENCIAL DISCRETA')

AGUSTÍN SOTO OTÁLORA

CONTROL DIGITAL

30

FIGURA 1.8

a k , x (k)    0,

para k  0, 1, 2, 3, ··· señal dis creta para k  0



X ( z )   a k z - k  sk  1  az 1  a 2 z 2  a 3 z 3    a k z  k k 0

Ec. (1)

Si multiplicamos la ecuación (1) en ambos miembros por az 1

az 1 sk  az 1  a 2 z 2  a 3 z 3    a k z  k  a k 1 z  ( k 1)

Ec. (2)

Si restamos la ecuación (2) - (1)

az 1 sk  sk  az 1  a 2 z 2  a 3 z 3    a k z  k  a k 1 z  ( k 1)  1  az 1  a 2 z 2  a 3 z 3    a k z  k

1

sk (az  1)  1  a

k 1  ( k 1)

z



 1  a k 1 z  ( k 1) sk  az 1  1

UNIVERSIDAD SURCOLOMBIANA

CAPÍTULO I

31

 1  a k a z 1 z  k  1  a  a z 1 1 z   1  a  a z 1 (0) lim sk  lim   k  k  az 1  1 az 1  1 az 1  1 lim sk  k 

X ( z) 

 1  (0) 1 1 z z     az 1  1 az 1  1 az  1 a  z z  a

1 z  1  az 1 z  a

La anterior relación se obtiene por Matlab de dos formas: Utilizando la definición de la transformada Z >> syms z k a

%se crean las variables simbólicas z, k, a

>> f=symsum(z^(-k)*a^(k),k,0,inf); >> pretty(f)

%se arregla la presentación z - ----a - z

Utilizando el comando ztrans

>> syms z k a

%se crean las variables simbólicas z, k, a

>> g=a^(k);

%se crea la función a transformar

>> f=ztrans(g);

%se calcula la transformada z

>> pretty(f)

%se arregla la presentación z -----------a (-1 + z/a) AGUSTÍN SOTO OTÁLORA

CONTROL DIGITAL

1.3.4

32

EXPONENCIAL: e -akT

(a = constante)

%GENERACION DE LA FUNCION EXPONENCIAL k = linspace (1,5,20);

x(k) = e

-2k

% define valores de k con % espaciamiento lineal

x = exp(-2* k);

% función exponencial

grid

% rejilla para gráfica

plot(k, x,'bo')

% gráfica x en función de k

xlabel('k');

% rotulo para eje x

ylabel('x(k)');

% rotulo para eje y

title('EXPONENCIAL DISCRETA')

FIGURA 1.9

UNIVERSIDAD SURCOLOMBIANA

CAPÍTULO I

33

para k  0, 1, 2, 3, ··· señal dis creta para k  0

e  akT , x (k)    0, 

X ( z )   e akT z - k  k 0

sk  1  e  aT z 1  e 2 aT z 2  e 3aT z 3    e  akT z  k Si multiplicamos la ecuación (1) en ambos miembros por

Ec. (1)

e  aT z 1

e  aT z 1sk  e  aT z 1  e 2 aT z 2  e 3aT z 3    e  aT ( k 1) z  ( k 1) Ec. (2) Si restamos la ecuación (2) - (1)

e  aT z 1 s k  s k  e  aT z 1  e 2 aT z 2  e 3aT z 3    e 2 akT z  ( k 1)  1  e  aT z 1  e 2 aT z 2  e 3aT z 3    e  akT z  k

sk (e

 aT

1

z  1)  e

 aT ( K 1)  ( K 1)

z

1



e aT ( K 1) z 8 K 1)  1 sk  e aT z 1  1

e aT ( K 1) z  ( K 1)  1 eaT 1 e z 1z   1 eaT (0) z 1 (0)  1 lim sk  lim   k  k  e aT z 1  1 e aT z 1  1 e aT z 1  1 lim s k  k 

X ( z) 

1 1 z   e  aT z 1  1 1  e  aT z 1 z  e  aT 1

1  e  aT z 1



z z  e  aT

La anterior relación se obtiene por Matlab de dos formas: Utilizando la definición de la transformada Z

AGUSTÍN SOTO OTÁLORA

CONTROL DIGITAL

>> syms z k a T

34

%se crean las variables simbólicas z k a T

>> f=symsum(z^(-k)*exp(-a*k*T),k,0,inf); >> pretty(f)

%se arregla la presentación z exp(a T) ---------------1 + z exp(a T)

Utilizando el comando ztrans >> syms z k a T

%se crean las variables simbólicas z k a T

>> g=exp(-a*k*T); %se crea la función a transformar >> f=ztrans(g);

%se calcula la transformada z

>> pretty(f)

%se arregla la presentación z -------------------------/

z

\

exp(-a T) |-1 + ---------| \

1.3.5

exp(-a T)/

SENOIDAL : sen(wkT)

%GENERACION DE LA FUNCION SENO: x(k) = sen(wkT) k = linspace(1,20);

% define valores de k con espaciamiento lineal

x = sin(k);

% función exponencial

grid

% rejilla para grafica UNIVERSIDAD SURCOLOMBIANA

CAPÍTULO I

35

plot(k, x,'bo')

% grafica x en función de k

xlabel('k');

% rotulo para eje x

ylabel('x(k) =seno(k)');

% rotulo para eje y

title('SENOIDAL DISCRETA')

FIGURA 1.10

sen(kT ) , x (k)   0, 

para k  0, 1, 2, 3, ··· señal dis creta para k  0

X(z) = Z[x(k)] = Z[sen(wkT)] Por al ecuación de Euler: Sen(wkT) = (1/2j) ( e jwkT - e – jwkT), reemplazando, AGUSTÍN SOTO OTÁLORA

CONTROL DIGITAL

36

Z[sen(wkT)] = Z[(1/2j)( e jwkT - e – jwkT)], aplicando la transf_z de la exponencial,

X ( z) 

1 1 1 1    2 j 1  e jT z 1 2 j 1  e  jT z 1

reemplazando las exponenciales :

z 1 sen(T ) X ( z)  1  2 z 1 cos(T )  z  2 1.4

PROPIEDADES Y TEOREMAS.

1.4.1 MULTIPLICACIÓN POR UNA CONSTANTE Z[a x(k)] = a Z[x(k)] = a X(z) 1.4.2 LINEALIDAD Si x(kT) = a f(kT) + b g(kT), entonces, X(z) = a F(z) + b G(z) 1.4.3

MULTIPLICACIÓN POR a k

Si y(kT) = a k x(kT), entonces, 

Z[a k y(kT)] =



a k x(kT) z – k

k 0





x(kT) (a

k 0

1.4.4

-1

z) – k = X(a

-1

z)

TEOREMA DE TRASLACIÓN

Si y(kT) = e

- akT

x(kT), entonces, 

Z[e

- akT

x(kT)] =

 k 0

e

-akT

x(kT) z – k

UNIVERSIDAD SURCOLOMBIANA

CAPÍTULO I

37





x(kT) (e

k 0

1.4.5

aT

z) – k = X(e

aT

z)

TEOREMA DEL CORRIMIENTO

Corrimiento hacia atrás: Z[x(k-n)T] = z – n Z[x(k)] = z

–n

X(z)

Corrimiento hacia adelante: n 1

Z[x(k+n)T] =

= z

n

z

X(z) - z

n

n

[ X(z) –

x(0) - z

n-1



x(kT)*z –k ]

k 0

x(1) - z

n-2

x(0) - z

2

x(2) - ········· - z x(n-1)

Ejemplo: Z[x(k+3)T] = z 1.4.6

3

X(z) - z

3

x(1) - z x(2)

SUMA DE FUNCIONES k

Sea y(k) =

 h0

x(h) , para k = 0,1, 2, ······

y(k) = x(0) + x(1) + x(2) + ········ + x(k-1) + x(k) y(k-1) = x(0) + x(1) + x(2) + ········ + x(k-1), restando estas dos expresiones, y(k) - y(k-1) = x(k), sacando Transf._Z, Y(z) – z

Y ( z) 

–1

Y(z) = X(z), entonces despejando Y(z) se tiene que:

1  X ( z) 1  z 1

AGUSTÍN SOTO OTÁLORA

CONTROL DIGITAL

.4.7

38

TEOREMA DEL VALOR INICIAL

Si el límite lim X(z) existe, entonces el valor inicial de x(k) = x(0) es igual a:

x(0)  lim x(k )  lim X ( z ) k 0

1.4.8

z 

TEOREMA DEL VALOR FINAL

El valor final de x(k), o sea, cuando k   (Si X(z) es estable) , es:



x()  lim x(k )  lim (1  z 1 ) X ( z) k 

z 1



Este se aplica si el sistema es estable. EJEMPLO 1-1 Encontrar la transformada Z de una función escalón de amplitud 4 y desfase en atraso de 3 periodos de muestreo.

Solución: x(kT) = 4*u(kT- 3T), asumiendo T = 1 por simplicidad, x(k) = 4u(k-3) Z[4u(k-3)] = 4Z[u(k-3)] = 4 z

-3

Z[u(k)]

aplicando teorema de corrimiento en atraso X(z) = 4 z

-3

(1/ (1-z

–1

)) = 4 / (z 3 – z 2)

EJEMPLO 1-2 Obtener la transformada Z de y(k) = 5

k–2

para k = 1, 2, 3, .... e igual a cero para k  0.

Solución: Sea x(k) = 5 k, entonces y(k) = x(k – 2) = 5 Z[y(k)] = Z[5 Z[5

k–2

k–2

] = Z[x(k -2)] = z

–2

k–2

Z[x(k)] = z

-2

Z[5 k ] = z – 2 * 1/(1 – 5 z

–1

)

] = 1 / (z 2 – 5 z)

EJEMPLO 1-3 Obtener la transformada Z de y(k) =k e

– 5k

para k = 1, 2, 3, .... e igual a cero para k ≤0.

Solución: UNIVERSIDAD SURCOLOMBIANA

CAPÍTULO I

39

Sea x(k) = k, entonces, X(z) = z

–1

/ (1 – z

–1

) 2 , y además, y(k) = e

– 5k

x(k ) ,

Aplicando teorema de traslación, Z[y(k)] = Z[k e

– 5k

] = X(e

5k

z )] , reemplazando 𝑧 = 𝑒 5𝑇 𝑧, en X(z) se tiene:

(𝑒 5𝑇 𝑧)−1 𝑒 −5𝑇 𝑧 −1 𝑌(𝑧) = = (1 − (𝑒 5𝑇 𝑧)−1 )2 (1 − 𝑒 −5𝑇 𝑧 −1 )2 EJEMPLO 1-4 Determinar el valor inicial x(0) de una señal si su transformada Z es igual a :

X ( z) 

(1  e 5k ) z 1 (1  z 1 )(1  e 5k z 1 )

Aplicando el Teorema de valor inicial,

(1  e 5k ) 1 x(0)  lim X ( z )  0 z  (1   1 )(1  e 5k  1 ) EJEMPLO 1-5 Determinar el valor final x(∞) de una señal si su transformada Z es igual a:

X ( z) 

1 1  (1  z 1 ) (1  e 5 z 1 )

 1  z 1  x()  lim x(k )  lim (1  z ) X ( z )  lim 1  1  5 1  k  z 1 z 1  1 e z 



1



EJEMPLO 1-6 Obtener la transformada Z de la figura dada. Tiempo de muestreo = 1.0

AGUSTÍN SOTO OTÁLORA

CONTROL DIGITAL

40

FIGURA 1.11

Si x(k) = (1/3)k (rampa de pendiente 1/3) y(k) = x(k) – x(k- 3), entonces, Y(z) = z[y(k)] = z[(1/3)k] –z

-3

z[(1/3)k]

1 z 1 1 z 3 * z 1 1 z 1 (1  z 3 ) Y ( z)       3 (1  z 1 ) 2 3 (1  z 1 ) 2 3 (1  z 1 ) 2 1.5 TRANSFORMADA Z INVERSA Con la transformada Z inversa se obtiene la señal discreta en los instantes de muestreo x(kT). Los siguientes son los métodos más utilizados para obtener la transformada Z inversa. 1.5.1

MÉTODO DE DIVISIÓN DIRECTA

El método consiste en arreglar la función X(z) en términos de z – 1 tanto el numerador como el denominador, dividir algebraicamente el numerador entre el denominador y su cociente mediante comparación con la definición de X(z) encontrar la señal x(kT). 

 x(k ) z

k



x(0) + x(1) z – 1 + x(2) z –2 + x(3) z – 3 + ······ + x(k) z – k

k 0

EJEMPLO 1-7 Obtener la transformada Z inversa x(k) de: UNIVERSIDAD SURCOLOMBIANA

CAPÍTULO I

X ( z) 

41

5 z  10 5 z  10  2 ( z  0.8)( z  0.2) z  z  0.16

multiplicando numerador y denominador por z

–2

,

5 z 1  10 z 2 X ( z)  1  z 1  0.16 z  2 dividiendo numerador entre denominador, se tiene, X(z) = 5 z

–1

+ 15 z

–2

+ 14.2 z

–3

+ 11.8 z

–4

+ ···

Comparando con la definición de X(z), se obtiene que: X(0) = 0, x(1) = 5, x(2) =15, x(3) = 14.2, x(4) = 11.8, ··· Si se quiere obtener más muestras de la señal, es mejor aplicar Matlab de esta forma:

% EJEMPLO 1-6: DE TRANSFORMADA Z INVERSA x = [1 zeros(1,40)];

% para k = 40 muestras

num = [0 5 10];

% coeficientes del numerador

den = [1 -1 0.16];

% coeficientes del denominador

y = filter(num,den,x) % obtención de las 40 muestras k = 0:40; plot(k,y,'ro',k,y,'-') xlabel('k') ylabel('y(k)') Los 40 resultados obtenidos de x(0) hasta x(39) son: 0

5.0000 15.0000 14.2000 11.8000

9.5280

7.6400

6.1155

4.8931

3.9146

3.1317

2.5054

2.0043

1.6035

1.2828

1.0262

0.8210

0.6568

0.5254

0.4203

0.3363

0.2690

0.2152

0.1722

0.1377

0.1102

0.0882

0.0705

0.0564

0.0451

0.0361

0.0289

0.0231

0.0185

0.0148

0.0118

0.0095

0.0076

0.0061

0.0048

0.0039

AGUSTÍN SOTO OTÁLORA

CONTROL DIGITAL

42

Cuya representación gráfica es la siguiente:

FIGURA 1.12 1.5.2 MÉTODO DE FRACCIONES PARCIALES Consiste el método en expandir la función X(z) en fracciones parciales con el fin de que queden términos más simples y luego encontrar a cada fracción la transformada Z inversa. EJEMPLO 1-8 Obtener la transformada Z Inversa de:

5 z 3  26 z 2  44 z  29 X ( z)  z 3  6 z 2  11z  6 Para representar esta función en fracciones parciales usamos el comando Matlab residue que encuentra los valores del vector r, del vector p y del término independiente k según la siguiente expresión:

X ( z) 

r r1 r  2  n  k z  p1 z  p 2 z  pn

Usando Matlab: num = [5

26

den = [1

6

44 11

29]; 6];

[r, p, k] = residue(num,den) % Los resultados son: r = [-2 -5 3],

p = [-3 -2 -1],

k = 5,

UNIVERSIDAD SURCOLOMBIANA

CAPÍTULO I

43

por tanto, las fracciones parciales quedan :

2 5 3   5 z  3 z  2 z 1

X ( z) 

Multiplicando por z-1 :

 2 z 1  5 z 1 3z 1 X ( z)    5 1  3z 1 1  2 z 1 1  z 1 Si se supone que: v(k) = ak Z[a

k–1

] = Z[v(k -1)] = z

–1

Z[v(k)] = z

-1

Z[a k ] = z – 1 * 1/(1 – a z

–1

)

entonces , x(k) = -2(-3)

k–1

-5(-2)k – 1 +3(-1)k – 1 + 5δk

donde δk es el delta kronecker (delta de Dirac en control continuo) cuya transformada Z es igual a 1. 1.5.3

MÉTODO DE LOS RESIDUOS

El método plantea que: m

x(kT) =



(residuos de X(z) z k - 1 en el polo z = zi )

i 1

x(kT) = k1 + k2 +······· + km (a) Si X(z) z k – 1 tiene un polo simple en z = zi , entonces, el residuo es :



ki  lim ( z  zi ) * X ( z ) *z k 1 z  zi



(b) Si X(z) z k – 1 tiene un polo múltiple en z = zi de orden q, entonces, el residuo es:

AGUSTÍN SOTO OTÁLORA

CONTROL DIGITAL

44



1  q1 q k 1 ki  lim ( z  z ) * X ( z ) * z i (q  1)! zzi z q1 EJEMPLO 1-9

CASO POLO SIMPLE. X(z) =

z (z − 2)(z − 3)(z − 1)

Se halla X(z)z k−1 𝑋(𝑧)𝑧 𝑘−1 =

𝑧𝑧 𝑘−1 (𝑧 − 2)(𝑧 − 3)(𝑧 − 1)

𝑋(𝑧)𝑧 𝑘−1 =

𝑧𝑘 , 𝑝𝑎𝑟𝑎 𝑘 = 0, 1, 2, … (𝑧 − 2)(𝑧 − 3)(𝑧 − 1)

𝑋(𝑧)𝑧 𝑘−1 tiene tres polos simples en z=2, z=3, z=1 𝑋(𝑘) = 𝑘1 + 𝑘2 + 𝑘3 Calculo de 𝑘1

(𝑧 − 2)𝑧 𝑘 2𝑘 𝑘1 = lim = = −2𝑘 𝑧→2 (𝑧 − 2)(𝑧 − 3)(𝑧 − 1) (−1)(1) (𝑧 − 3)𝑧 𝑘 3𝑘 3𝑘 𝑘2 = lim = = 𝑧→3 (𝑧 − 2)(𝑧 − 3)(𝑧 − 1) (1)(2) 2 (𝑧 − 1)𝑧 𝑘 1𝑘 1 𝑘3 = lim = = 𝑧→1 (𝑧 − 2)(𝑧 − 3)(𝑧 − 1) (−1)(−2) 2

𝑋(𝑘) = −2𝑘 +

3𝑘 1 + 2 2

UNIVERSIDAD SURCOLOMBIANA



CAPÍTULO I

45

POR MATLAB

num=[0 0 1 0];% coeficientes del numerador hay que colocar los ceros para hallar Xz1 den=[1 -6 11 -6];% coeficientes del denominador Xz=tf(num,den,-1); % función de transferencia en términos de z pole(Xz)% obtención de polos % la función tiene

polos simple en z =2, z=3 z=1

% se debe reconstruir la función para aplicar el % comando limit de Matlab syms z Xz1 = sum(num.*[z^3 z^2 z 1])/sum(den.*[z^3 z^2 z 1]) % a) cálculo del residuo para k1,k2,k3 plos simples syms z k k1 =limit ((z-2)*Xz1*z ^(k-1),z,2) k2=limit((z-3)*Xz1*z^(k-1),z,3) k3=limit((z-1)*Xz1*z^(k-1),z,1) %X(k)=k1+k2+k3

CASO POLO MULTIPLE.

3z 2  9 z X ( z)  3 z  1.8 z 2  1.05 z  0.2 X(k)=K1+K2 Hallamos K1 que corresponde al polo simple

 3z 2  9 z k 1  k1  lim z  0.8* * z  2 z 0.8    z  0 . 8 z  0 . 5   AGUSTÍN SOTO OTÁLORA

CONTROL DIGITAL

46

 3z  9 k k1  lim  * z  z 0.8  z  0.52   3 * 0.8  3 * 0.8 3 *  2.2 * 0.8  6.6 * 0.8 k1    0.09 0.8  0.52 0.32 k

k

k

k1  73.30.8

k

Ahora hallamos K2 que corresponde al polo múltiple



1  q1 q k 1 ki  lim ( z  z ) * X ( z ) * z i (q  1)! zzi z q1

𝑘2 =



1 𝜕 (𝑧 − 0.5)2 (3𝑧 − 9)𝑧 𝑘 lim [ ] (2 − 1)! 𝑧→0.5 𝜕𝑧 (𝑧 − 0.5)2 (𝑧 − 0.8)

𝑑 3𝑧 𝑘+1 − 9𝑧 𝑘 𝑘2 = lim ( ) 𝑧→2 𝑑𝑧 (𝑧 − 0.8) (𝑧 − 0.8)[3(𝑘 + 1)𝑧 𝑘 − 9𝑘𝑧 𝑘−1 ] − [3𝑧 𝑘+1 − 9𝑧 𝑘 ](1) 𝑘2 = (𝑧 − 0.8)2 (𝑧 − 0.8)[3𝑘𝑧 𝑘 + 3𝑧 𝑘 − 9𝑘𝑧 𝑘−1 ] − 3𝑧 𝑘+1 + 9𝑧 𝑘 𝑘2 = (𝑧 − 0.8)2 (−0.3)[3𝑘(0.5)𝑘 + 3(0.5)𝑘 − 9𝑘(0.5)𝑘−1 ] − 3(0.5)𝑘+1 + 9(0.5)𝑘 𝑘2 = (−0.3)2 2.7𝑘(0.5)𝑘 −0.9𝑘(0.5) − 0.9(0.5) + − 3(0.5)(0.5)𝑘 + 9(0.5)𝑘 0.5 𝑘2 = 0.09 𝑘

𝑘

UNIVERSIDAD SURCOLOMBIANA

CAPÍTULO I

47

−0.9𝑘(0.5)𝑘 − 0.9(0.5)𝑘 + 5.4𝑘(0.5)𝑘 − 1.5(0.5)𝑘 + 9(0.5)𝑘 𝑘2 = 0.09 𝑘2 = −10𝑘(0.5)𝑘 − 10(0.5)𝑘 + 60𝑘(0.5)𝑘 − 16.66(0.5)𝑘 + 100(0.5)𝑘 𝑘2 = 50𝑘(0.5)𝑘 + 73.33(0.5)𝑘 Ahora resolvamos el mismo ejercicio por MAT-LAB % EJEMPLO 1-9 : PROGRAMA EN MATLAB num=[0 3 -9 0];

% coeficientes del numerador

den=[1 -1.8 1.05 -0.2];% coeficientes del denominador Xz=tf(num,den,-1);

% función de transferencia en % términos de z

pole(Xz)

% obtención de polos

% la función tiene un polo simple en z = 0.8 y un polo % doble en z = 0.5 % se debe reconstruir la función para aplicar el % comando limit de Matlab syms z Xz1 = sum(num.*[z^3 z^2 z 1])/sum(den.*[z^3 z^2 z 1]);

% a) cálculo del residuo para el polo simple syms z k k1 =limit ((z-0.8)*Xz1*z ^(k-1),z,0.8) % k1 = -220/3*4^ k/(5^ k) AGUSTÍN SOTO OTÁLORA

CONTROL DIGITAL

48

% b) cálculo del residuo para el polo doble h=diff((z-0.5)^2*Xz1*z^(k-1),1)

% primera derivada

k2=limit(h,z,0.5) %K2=1/3*(220+150*k)/(2^k)

x(k) = k1 + k2 =

1.5.4

-220/3*4^ k/(5^ k) + 1/3*(220+150*k)/(2^k)

UTILIZANDO MATLAB PARA HALLAR Z INVERSA

Para calcular una transformada Z inversa se utiliza el comando iztrans >> syms z k

%se crean las variables simbólicas z,k

>> g=z/(z-1);

%se crea la función a transformar

>> f=iztrans(g);%se calcula la transformada z inversa >> pretty(f)

%se arregla la presentación

>> syms z k

%se crean las variables simbólicas z,k

>> g=z/(z-1)^2;

%se crea la función a transformar

>> f=iztrans(g); %se calcula la transformada z inversa >> pretty(f)

%se arregla la present

>> syms z k a

%se crean las variables simbólicas z,k

>> g=z/(z-a);

%se crea la función a transformar

>> f=iztrans(g);

%se calcula la transformada z

>> pretty(f)

%se arregla la presentación

1.6 ECUACIONES EN DIFERENCIA Considérese un sistema discreto LTI (Lineal e invariante en el tiempo) dado por la ecuación en diferencias: x(k) + a1 x(k-1) + a2 x(k-2) +··· +an x(k-n) = b0 u(k) + b1 u(k-1) + b2 u(k-2) +····+bn u(k-n)

UNIVERSIDAD SURCOLOMBIANA

CAPÍTULO I

49

donde u(k) es la entrada al sistema y x(k) es la salida. La forma de solucionar esta ecuación en diferencia consiste en calcular la transformada Z, luego aplicar las condiciones iniciales dadas y por último obtener la transformada Z inversa. Se debe recordar que: Z[x(k-n)T] = z – n Z[x(k)] = z Z[x(k+n)T] = z

n

–n

n 1



[ X(z) –

X(z),

y,

x(kT)*z –k ] =

k 0

z

X(z) - z

n

n

x(0) - z

n-1

x(1) - z

n-2

x(2) - ········· - z x(n-1)

EJEMPLO 1-10 Resolver la siguiente ecuación en diferencias. x( k +2) + 5x(k +1) + 6x(k) = 0 Condiciones iniciales: x(0) = 0, x(1) = 1

Solución : a) Aplicando Transf._Z, se tiene: [z

2

X(z) - z

2

x(0) - z x(1)] + 5[z X(z) - z x(0)] + 6 X(z) = 0

b) Sustituyendo condiciones iniciales: [z

2

X(z) - z =z

2

2

(0) - z (1)] + 5[z X(z) - z (0)] + 6 X(z) = 0

X(z) - z + 5z X(z) + 6 X(z) = 0

despejando X(z):

X ( z) 

z z z   z  5z  6 z  2 z  3

(fracciones parciales)

2

c) Obtener la transf_z inversa : Sabiendo que, z[ ak ] = 1/ (1 – a z

X ( z) 

–1

), entonces ,

1 1  1 1  2z 1  3z 1

Por tanto su Transf._Z inversa es: x(k) = (-2)k – (-3)k

AGUSTÍN SOTO OTÁLORA

CONTROL DIGITAL

50

EJEMPLO 1-11 Resolver la siguiente ecuación en diferencias. a) Aplicando Transf._Z, se tiene: x(k) + 5x(k-1) + 6x(k-2) = u(k), unitario

donde u(k) es el escalón

Solución : X(z) + 5 z

–1

X(z) + 6 z

-2

X(z) = U(z), pero U(z) = 1 / (1 – z- 1)

Despejando X(z) :

1 z3 X ( z)   (1  z 1 )(1  5 z 1  6 z  2 ) ( z  1)( z 2  5 z  6) b) Obtener la transf_z inversa: % Aplicando Matlab para encontrar las fracciones parciales Xz = zpk([0

0

0

], [1

-2

-3], 1, -1);

Xz = tf(Xz) pole(Xz) % tiene polos simples en : -3.0000

-2.0000

1.0000

[num,den]=tfdata(Xz,'v') [r,p,k]=residue(num,den) % r =

-6.7500

2.6667

0.0833

% p =

-3.0000

-2.0000

1.0000

% k =

1

Con base en lo anterior,

X ( z) 

z  6.75 z  2.6667 z  0.0833   1 z 3 z2 z 1

Multiplicando por z-1 :

X ( z) 

1  6.75 z 1 1  2.6667 z 1 1  0.0833z 1   1 1  3z 1 1  2 z 1 1  z 1

Sabiendo que, z[ak] = 1/ (1 – a z

–1

), entonces ,

x(k) = (-3)k + 6.75(-3)k - 1 + (-2)k - 2.6667(-2)k – 1 + (1)k – 0.0833(1)k – 1 + δk UNIVERSIDAD SURCOLOMBIANA

CAPÍTULO I

51

= (-3)k + 6.75(-1/3)(-3)k + (-2)k – 2.6667(-1/2) (-2)k + 1 - 0.0833 + δk x(k) = -5/4*(-3)k + 7/3*(-2)k + 0.9167 + δk

EJEMPLO 1-12 Resolver la siguiente ecuación en diferencias. x (k + 2) + 0.5x (k + 1) + 0.2x(k) = u(k + 1) + 0.3u(k),

(1)

Condiciones iniciales: x(k) = 0, para k  0 y u(k) = 0, para k  0 y además, u(0) = 1.5, u(1) = 0.5, u(2) = -0.5, u(k) =0 para k = 3, 4, 5, ·······

Solución : a) Aplicando Transf._Z, se tiene: [ z 2 X(z) – z2 x(0) – z x(1) ] + 0.5[ z X(z) – z x(0)] + 0.2 X(z) = [z U(z) –z u(0)] + 0.3 U(z) (2) b) Sustituyendo condiciones iniciales: Como no se conoce x(1), se debe encontrar reemplazando k por -1 en la ecuación (1), x(-1+2) + 0.5x(-1+1) + 0.2x(-1) = u(-1+1) + 0.3u(-1),

entonces,

x(1) + 0.5 x(0) + 0.2 x(-1) = u(0) + 0.3 u(-1), de las condiciones iniciales, x(0)=0, x(-1)= 0, u(0) = 1.5, u(-1) = 0, reemplazando se tiene que, x(1) = 1.5 encontrado x(1) se sustituyen condiciones iniciales en (2), [ z 2 X(z) – z (1.5) ] + 0.5[ z X(z) ] + 0.2 X(z) = [z U(z) - z (1.5)] + 0.3 U(z) (2) Despejando X(z):

X ( z) 

z  0.3 U ( z) z 2  0.5 z  0.2 

U ( z )   u (k ) z  k u (0)  u (1) z 1  u (2) z  2  u (3) z 3    k 0

k=0 = 1.5 + 0.5 z – 1 - 0.5 z – 2 , reemplazando, AGUSTÍN SOTO OTÁLORA

CONTROL DIGITAL

52

( z  0.3)(1.5 z 2  0.5 z  0.5) 1.5 z 3  0.95 z 2  0.35 z  0.15 X ( z)   ( z 2  0.5 z  0.2) z 2 z 4  0.5 z 3  0.2 z 2 Para obtener la señal discreta por Matlab:

% OBTENER DIEZ VALORES DE LA SEÑAL num = [0 1.5 0.95 -0.35 -0.15]; den = [1 0.5 0.2 0 0]; u = [1 zeros(1,10)]; xk = filter(num,den,u) k = 0:10; stem(xk, k)

% grafica la señal muestreada

Valores x(k) dados por Matlab: x(0)= 0, x(1)= 1.5000, x(2)= 0.2000, x(3)=-0.7500, x(4)= 0.1850, x(5)= 0.0575 x(6)= -0.0658, x(7)= 0.0214, x(8)= 0.0025, x(9)= -0.0055, x(10)= 0.0023

Ejemplo 1-13 canibales

y[k  2]  1.5 y[k  1]  0.8 y[k ]  u[k  2]

(1)

Considerando la transformada Z de una secuencia desplazada en el tiempo, la TZ de (1) será,

Y ( z ) z 2  1.5Y ( z ) z  0.8Y ( z )  U ( z ) z 2 De esta manera, se tiene que, 𝑌(𝑧)[𝑧 2 − 1.5𝑧 + 0.8] = 𝑈(𝑧)𝑧 2 De donde, 𝑌(𝑧) =

𝑈(𝑧)𝑧 2 [𝑧 2 − 1.5𝑧 + 0.8]

(2)

Ahora bien, definiendo la transformada Z de la función de excitación (entrada),

UNIVERSIDAD SURCOLOMBIANA

CAPÍTULO I

𝑈(𝑧) =

2𝑧 𝑧−1

53

(3)

Considerando (3) en (2), 𝑌(𝑧) =

𝑧2 2𝑧 ∙ 2 [𝑧 − 1.5𝑧 + 0.8] 𝑧 − 1

(4)

Ahora bien, para que la función racional (4) sea una fracción propia, y se puede realizar su descomposición en fracciones parciales, se transformara (4) como sigue:

𝑌(𝑧) 𝑧2 2 = 2 ∙ [𝑧 𝑧 − 1.5𝑧 + 0.8] 𝑧 − 1

(5)

𝑧2 2 𝐴 𝐵 𝐶 ∙ = + + 2 [𝑧 − 1.5𝑧 + 0.8] 𝑧 − 1 𝑧 − 1 𝑧 − 𝛼1 𝑧 − 𝛼2

(6)

Donde 𝛼1 y 𝛼2 serán numero complejos conjugados, considerando el orden del polinomio de z. Realizando la descomposición a través de MatLab se tiene, >> num = [2 0 0 ]; >> den = [1 -2.5 2.3 -0.8]; >> [r ,p, k] = residue(num,den); r = 6.6667 -2.3333 - 1.8810i -2.3333 + 1.8810i p = 1.0000 0.7500 + 0.4873i AGUSTÍN SOTO OTÁLORA

CONTROL DIGITAL

54

0.7500 - 0.4873i k = [] Donde la estructura [r, p, k] corresponde a los coeficientes –r-, los polos de la fracción racional –p-, y los términos enteros –k-.De esta manera se tiene que la expresión (5) se transforma en:

𝑌(𝑧) 6.6667 −2.3333 − 𝑗1.8810 −2.3333 + 𝑗1.8810 = + + 𝑧 𝑧−1 𝑧 − (0.75 + 𝑗0.4873) 𝑧 − (0.75 − 𝑗0.4873)

Determinando la transformada inversa Z de la expresión

(7)

(7) se tiene:

𝑌[𝑛] = 6.6667 + (−2.3333 − 𝑗1.8810)(0.75 + 𝑗0.4873)𝑛 + (−2.3333 + 𝑗1.8810)(0.75 – 𝑗0.4873)𝑛 M-FILE asociado function y = tz(x) num = [2 0 0 ]; den = [1 -2.5 2.3 -0.8]; [r ,p, k] = residue(num,den);

y = (r(1))*(p(1).^x) + (r(2))* (p(2).^x) + (r(3))*(p(3).^x); end

EJECUCION EN EL WORKSPACE

>> x = 0:1:15; >> tz(x)

UNIVERSIDAD SURCOLOMBIANA

CAPÍTULO I

55

FIGURA 1.13 TRANSFORMADA Z MODIFICADA. Los comportamientos entre los puntos de muestreo pueden ser investigados utilizando la transformada Z modificada. Esta es la transformada Z ordinaria, solamente retrasada mT segundos, la cual es una fracción del periodo de muestreo, ya que 0