Modelado de Sistemas de Control Digital

DR. FERNANDO ORNELAS TELLEZ

25 de octubre de 2016

Resumen
El alumno comprenderá los métodos de análisis y diseño de los sistemas de control digital, así como la
realización de algoritmos de control.

1. Introducción

Dado que el sistema de control digital sólo toma datos en forma discreta del sistema continuo, entonces la
computadora considera al sistema continuo como si éste fuese un sistema discreto, por lo tanto surge la necesidad
de modelarlo. Es evidente que al muestrear el sistema continuo, sólo se analiza y/o procesa la información discreta
disponible en el instante de muestreo, ignorando lo que pueda ocurrir entre dichos instantes, haciendo inevitable
la perdida de información. Tanto en la teoría como en la práctica, se considera que la información perdida no será
crucial para los fines de control, en el entendido que se está satisfaciendo cabalmente el Teorema de Muestreo
(o se muestrea mucho mas rápido que la frecuencia de Nyquist), o que el error por el proceso de muestreo entre
el sistema continuo y el discreto, es lo suficientemente pequeño para poder diseñar controladores basados en las
representaciones discretas.

1.1. Discretizacion de Controladores Lineales Continuos

Consideraciones sobre la versión discreta del esquema de control:

Se tiene perdida de información debido a que solo se tiene información disponible en los instantes de
muestreo;

1
 Perdida de información debida a la cuantización (ciertos valores son redondeados a una cantidad que se
pueda representar en un numero de bits específico);
Retraso en la respuesta del controlador;
Dependiendo del tipo de discretizacion;

Velocidad de muestreo limitados.

Aunque se tienen varias desventajas como las anteriores, su implementación es fácilmente realizable en cualquier
computador digital apropiado; Se tiene la flexibilidad de que un rediseño en el esquema de control, no implica
remplazar ningún componente, solo basta re-programar; pueden resultar en prototipos de control económicos.
Una forma simple de realizar aproximado digital de un controlador a partir de un controlador analógico es
usar el método de Euler, considerando
x(k + 1) x(k)
ẋ(k) ⇠
=
T
donde T = tk+1 tk es el tiempo de muestreo, tk es el tiempo en el instante k, x(k) es el calor de x en el tiempo
tk . Esta aproximación1 se utiliza para pasar de una ecuación diferencial a un conjunto de ecuaciones que puedan
ser resueltas en una computador digital. Esas ecuaciones se les llama ecuaciones en diferencia. Ha de tenerse
cuidado dicha aproximación, pues debe considerarse la rapidez de la dinámica o rapidez del sistema continuo
y así fijar el periodo de muestreo T . Por ejemplo, un sistema que tiene un ancho de banda del orden de los
100 Hz, el periodo de muestreo puede ser del orden de los 10 msec y los errores de la aproximación resultaran
bastante pequeños para fines de control.

1.1.1. Ecuaciones en Diferencias

Considere que a k como una variable independiente representando el tiempo discreto (segundos, horas, días,
etc.) que ha transcurrido desde un momento inicial k = 0. Así, {y0 , y1 , ... , yk , ...} es una sucesión donde yk
corresponde a un valor de y en el instante k.
Definición 1. Dada una sucesión {xk } cuyos primeros términos son x0 , x1 , x2 , ..., presentamos como Ecua-
ción en Diferencias a toda ecuación que relaciona términos de esta sucesión.
Definición 2. Llamamos ecuación en diferencias a una expresión del tipo

F (yk+n , yk+n 1, ... , yk+1 , yk ) = 0.

Usa solución de la misma, es una sucesión de calores de y que cumpla tal expresión.
El orden de una ecuación en diferencias la define la diferencia entre el mayor y el menor de los índices de la
ecuación. Por ejemplo el orden de 2xk+3 + 3xk = 0, es 3.
El conjunto de todas las soluciones recibe el nombre de solución general, misma que estará en función de
parámetros por definir a partir de las condiciones iniciales, obteniendose soluciones particulares.

1.1.2. Ejemplo: Ecuación en Diferencias Utilizando el Método de Euler

Obtenga un controlador discreto a partir de uno digital cuya función de transferencia es [4]:

U (s) s+a
D(s) = =K .
E(s) s+b

Solución: La versión discreta es:

u(k + 1) = (1 b T ) u(k) + K (a T 1) e(k) + K e(k + 1).

1A esta aproximación a menudo se le encuentra con el nombre de regla rectangular hacia adelante.

2
 Tarea: Determine un controlador digital utilizando la aproximación de Euler para el compensador de ade-
lanto de fase [4]:
s+2
D(s) = 70
s + 10
para la planta o sistema
1
G(s) =
s (s + 1)
usando una frecuencia de muestreo de 20 Hz y otra de 40 Hz. Realice su simulación y compare tanto la respuesta
del sistema en el tiempo del controlador continuo como el discreto obtenido mediante la aproximación.

1.2. Selección del Periodo de Muestreo

La selección del periodo de muestreo es un aspecto muy importante en sistemas controlados por un compu-
tador digital. Un periodo de muestreo largo podría hacer imposible la reconstrucción de la señal en tiempo
continuo, mientras que un periodo de muestreo muy corto, incrementará la carga computacional en el dispositi-
vo digital utilizado [2]. En esta sección se aborda la selección del periodo de muestreo en función de la respuesta
de un sistema de control en tiempo continuo a lazo abierto. Dado lo anterior, se considerará el tiempo de subida
de la salida de un sistema como marco de referencia para establecer un adecuado periodo de muestreo [2].
Introduciendo Nr como el número de muestras por tiempo de subida, entonces
Tr
Nr =
T
donde Tr es el tiempo de subida de un sistema. Para sistemas de primer orden, el tiempo de subida Tr se
establece igual a la constante de tiempo ⌧ . Es entonces razonable elegir Nr entre 4 y 10, esto es Nr = 4 ⇠ 10,
⌧ ⌧
por lo tanto el periodo de muestreo para un sistema de primer orden es T = = . Para un sistema de
Nr 4 ⇠ 10
segundo orden (el cual se considera que tiene tiene polos con parte imaginaria) con factor de amortiguamiento
0 < ⇣ < 1 y frecuencia natural !n en rad/s, el tiempo de subida se define como

Tr = !n 1 e / tan

donde ⇣ = cos , o bien, para un valor de ⇣ alrededor de ⇣ ⇡ 0.701, de forma aproximada se tiene que Tr t
0.2 ⇠ 0.6
[2]. También, de manera aproximada [4], se tiene que
!n
1.8
Tr t
!n
Tr
y por tanto se podría elegir un periodo de muestreo adecuado al seleccionarlo como Ts = . Otra aproximación
Nr
para el tiempo de subida es el encontrado en [5], con valor
1 0.4167⇣ + 2.917⇣ 2
Tr = , 0 < ⇣ < 1.
!n

Otros métodos heurísticos se consideran en [9] para la determinación del periodo de muestreo apropiado,
así, se considera como regla practica el muestrear de ocho a diez veces durante el tiempo de levantamiento de
la respuesta de un sistema sobre-amortiguado ante una entrada tipo escalón, mientras que para un sistema
subamortiguado (oscilatorio), se debe considerar muestrear de 8 a 10 veces durante un ciclo de las oscilaciones
senoidales amortiguadas a la salida del sistema de control. Para evaluar la ubicación de los polos del sistema
muestreado a partir del sistema continuo, se puede
p considerar la transformación z = eT s , por ejemplo. para un
sistema que tiene un polo en s = ⇣!n + j!n 1 ⇣ 2 , el polo correspondiente discreto estará en
⇣ p ⌘
T ⇣!n +j!n 1 ⇣ 2
z = eT s = e .

De la relación anterior podemos destacar que la selección adecuada de T será de suma importancia en la
estabilidad del sistema, ya que para valores grandes de T , los polos en el plano complejo z pueden estar fuera

3
del circulo unitario y por tanto ser inestable el sistema. Para un factor de amortiguamiento alrededor de ⇣ ⇡ 0.7,
se obtiene que el periodo de muestreo se puede calcular como !n T ⇡ 0.2 0.6 [2].
Tarea: Determinar la versión discreta del PID e implementarla como un controlador digital. Determine un
periodo de muestreo apropiado y realice su simulación. Para este ejemplo, considere la función de transferencia
de un servomotor [4]:
360000
G(s) =
(s + 600) (s + 600)
y una función de transferencia para el PID continuo como
1
D(s) = K(1 + + TD s)
TI s
donde K = 3.2, TD = 0.0011 segundos y TI = 0.003 segundos.

2. Obtención de Modelos Discretos a partir de Modelos Continuos

2.1. Ecuación de Estado y su Solución
La solución de la ecuación de estado no homogénea (sistema lineal invariante en el tiempo):

ẋ = Ax + Bu, x(0) = x0 (1)

y = Cx+Du (2)
puede obtenerse mediante el enfoque de la transformada de Laplace. La transformada de Laplace de esta última
ecuación produce:
s X(s) x(0) = A X(s) + B U (s)
o
(s I A) X(s) = x(t0 ) + B U (s)
1
Pre-multiplicando ambos miembros de esta última ecuación por (s I A) , obtenemos
1 1
X(s) = (s I A) x(0) + (s I A) B U (s)
⇥ ⇤ ⇥ ⇤
= L eAt x(0) + L eAt B U (s)

La transformada inversa de Laplace se obtiene mediante la integral de convolución, donde L 1

{F (s) G(s)} =
Rt
(f ? g) (t) = 0 f (t ⌧ ) g (⌧ ) d⌧ , de modo siguiente:
Z t
At
x(t) = e x(0) + eA(t ⌧)
B u(⌧ ) d⌧ (3)
0
En caso de que el tiempo inicial sea diferente de cero, entonces se considera t0 en lugar de 0, y la solución
anterior se modifica a
Z t
x(t) = eA(t t0 )
x(t0 ) + eA(t ⌧)
B u(⌧ ) d⌧ (4)
t0
⇥ ⇤ 1
Notar que se utilizo la definición L eAt = (s I A) . En ocasiones la ecuación anterior se expresa en terminaos
de la matriz de transición: (t, t0 ) = eA(t t0 ) . Note que la ecuación anterior no se puede resolver a menos que
se hagan ciertas hipótesis sobre la entrada u(t).

4
2.1.1. Solución de la Ecuación de Estado en el Dominio del Tiempo

Considere el sistema descrito por (16) y (2).

Recordando la siguiente relación
d At
e = AeAt = eAt A.
dt
Pre-multiplicando (16) por e At
se tiene
At At At
e ẋ(t) = e A x(t) + e B u(t)

o bien
At At At
e ẋ(t) e A x(t) = e B u(t)
Entonces
d At At
e x(t) = e B u(t)
dt
Integrando la expresión anterior se obtiene
t
Z t
At A⌧
e x(t) = e B u(⌧ ) d⌧
t0 t0

y por tanto
Z t
At t0 A⌧
e x(t) e x(t0 ) = e B u(⌧ ) d⌧
t0

Considerando que la inversa de e At es eAt , y que e0 = I (la matriz identidad de dimensión apropiada), la
solución de la ecuación de estado (16) es la ecuación (4). De (2), es inmediato expresar la salida como
Z t
A(t t0 )
y(t) = Ce x(0) + C eA(t ⌧)
B u(⌧ ) d⌧ + Du(t).
t0

En matlab es posible obtener la versión discreta a partir de una continua en representación de espacio de
estados, considerando un periodo de muestreo T s. La función es c2d().

2.1.2. Discretización Exacta de un Sistema LTI Continuo

Considerando un periodo de muestreo constante T , donde se tiene que tk = k T , el sistema discretizado

resulta en
x((k + 1)T ) = x(kT ) + u(kT )
y(kT ) = C x(k T ) + D u(k T )
donde
= eA T
Z T
= eA⌧ d⌧ B
0

bajo la consideración de que u es constante entre instantes de muestreo (por ejemplo u como salida de un
retenedor de orden cero). El modelo anterior se le conoce como muestreo basado en el retenedor de orden cero
[2]. Un análisis semejante esta en [10] . Otro análisis parecido esta en [7].
En [9] se da detalladamente la discretización de sistemas continuos representados en espacio de estados, así
como diferentes ejemplos .

5
2.2. Solución del sistema lineal de primer orden
Considere el sistema lineal invariante en el tiempo con entrada u(t) y salida y(t) cuya función de transferencia
es
Y (s) b
G(s) = = (5)
U (s) s+a
donde Y (s) y U (s) es la transformada de Laplace de u(t) y y(t), respectivamente, con a y b constantes reales.
La ecuación diferencial respectiva al sistema (5) resulta en

ẏ(t) + ay(t) = bu(t) (6)

De acuerdo a la solución dada en (4), que corresponde a un sistema lineal invariante en el tiempo, tenemos
que la solución para (6) es
Z t
y(t) = e a(t t0 ) y(t0 ) + e a(t ⌧ ) b u(⌧ ) d⌧ (7)
t0

2.3. Modelo discreto exacto a partir de un sistema lineal con entrada constante
Considere que para la solución dada en (7), se tiene una entrada tipo escalón unitario u(t) = 1, una condición
inicial y(t0 ) = y(kT ), y remplazando t por (k + 1) T , entonces se obtendrá la solución y(t) en tiempo discreto
como y((k + 1)T ), en función de y(k T ):
Z (k+1)T
aT a(kT +T ⌧)
y((k + 1)T ) = e y(k T ) + e b u(⌧ ) d⌧
kT
✓ ◆ Z (k+1)T
aT 1 a(kT +T ⌧)
= e y(k T ) + bu(k T ) e a d⌧
a kT

que finalmente se expresa como

aT b aT
y((k + 1)T ) = e| {z } y(k T ) + a 1 e u(kT )
↵:= | {z }
:=
= ↵ y(k T ) + u(k T )

Ha de resaltarse que la solución obtenida es exacta en cada instante de muestreo, es decir, la solución
del modelo continuo y el discreto van a coincidir solo en los instantes en los cuales se realiza el muestreo,
considerando que el sistema es lineal, con entrada constante (o constante en cada instante de muestreo como
resultado del retenedor de orden cero).

2.3.1. Esquemas para simulación

Implementación en Simulink de y((k + 1)T ) = ↵ y(kT ) + u(kT ).

6
2.3.2. Ejercicio

Considere el sistema lineal (5) con a = 1, b = 1 y T = 0.5. Obtenga su modelo discreto por retenedor de
orden cero e implemente su simulación. Compare la solución del sistema continuo y(t) y del sistema discreto
y(kT ). De igual forma, compare la solución con la versión aproximada de Euler.

2.4. Propiedades Importantes de Sistemas Discretos

2.4.1. Teorema de Traslación (corrimiento) Real

Corrimiento a la derecha (retraso en el tiempo o corrimiento hacia atrás): es denotado por q 1

óz 1
así, un
corrimiento unitario se puede escribir como
1
z f (kT ) = f ((k 1)T ).

Deforma general se tiene:

n
z f (kT ) = f ((k n)T ).
Mientras que el operador de corrimiento hacia adelante u operador de adelanto z

z f (kT ) = f ((k + 1)T ).

Deforma general se tiene:

z n f (kT ) = f ((k + n)T ).

2.5. Aproximación de la derivada por diferencias finitas

Considere el problema de aproximar la derivada de una función x(t) conociendo únicamente datos de la señal
en instantes de tiempo consecutivos, como se muestra en la siguiente Figura

Entonces una aproximación a la derivada en el tiempo de x(t) se puede obtener como

d x(t)
ẋ =
dt
x
⇡
t
x(t + t) x(t)
=
t
x(kT + T ) x(kT )
=
T
A ésta aproximación de la derivada se le nombra aproximación por diferencias hacia adelante, misma que en
términos del operador de corrimiento se puede expresar de la siguiente manera:
q 1
ẋ(kT ) ⇡ x(kT ).
T

7
 De manera semejante, se puede definir la aproximación por diferencias hacia atrás como
d x(t)
ẋ =
dt
x
⇡
t
x(t) x(t t)
=
t
x(kT ) x(kT T )
T
mismo que en términos del factor de corrimiento se expresa como
1
1 q
ẋ(kT ) ⇡ x(kT ).
T

No obstante, la aproximación por diferencias hacia atrás presenta un error de aproximación mayor con
respecto al de diferencias hacia adelante.

2.6. Aproximación de la integral (Regla trapezoidal o de Tustin)

Planteando un sistema dinámico mediante integrales en lugar de ecuaciones diferenciales, un aproximador de
la operación de integral es la regla trapezoidal o de Tustin, que tiene la idea básica de obtener una aproximación
del área bajo la curva entre dos instantes de muestreo.

Rt
Definiendo la integral bajo la curva como I(t) = 0 x(⌧ ) d⌧ , por tanto considerando un tiempo discreto
tk = kT y un tiempo tk+1 = (k + 1)T , la integral bajo la curva vendrá dada como

x [(k + 1)T ] + x [kT ]

I [(k + 1)T ] I [kT ] ⇡ T
2
misma que en términos del operador de corrimiento se puede escribir

T q+1
I(kT ) ⇡ x(kT )
2 q 1
De esta forma
Z t
1
X(s) = x(⌧ ) d⌧
s 0
= I(kT )

T q+1
⇡ x(kT )
2 q 1
Por lo tanto 
1 T q+1
⇡
s 2 q 1

8
mismo resultado que puede ser utilizado para aproximar la derivada considerando la operación de inversa a cada
lado de la aproximación, esto es 
2 q 1
s⇡ .
T q+1

Ejemplo: Realice las aproximaciones a la derivada anteriores (aproximación hacia adelante, hacia atrás y
la de Tustin) al sistema lineal de primer orden

ẋ + a x = b u, x(t0 ) = 0

y compare los resultados obtenidos considerando un periodo de muestreo T = 0.1. ¿Cuál es mejor? ¿Cuál es
mas complicado?

2.7. Modelos discretos de sistemas continuos de orden superior

Esta sección trata lo referente a la discretización aproximada de sistemas de orden superior (mayor a grado
uno), mismos que aparecen en sistemas típicos de control. La solución a esto se puede tratar considerando las
aproximaciones a la variable compleja s, la cual indica la operación de derivada en el dominio de la frecuencia.
Resumiendo las aproximaciones a la derivada, se puede plantear la siguiente tabla:

Método de Aproximación Aproximación

q 1
Por diferencia hacia adelante s!
T
1 q 1
Por diferencia hacia atrás s!
T
2 q 1
Tustin s!
T q+1

Ejemplo: Obtener la aproximación discreta del siguiente sistema de segundo orden

ẍ + a1 ẋ + a0 x = b u

utilizando cada uno de los aproximadores de derivada.

Sol. Aplicando la transformada de Laplace a la ecuación, se llega a que
⇥ 2 ⇤
s + a1 s + a0 Y (s) = b U (s)

por lo que usando la aproximación de diferencias hacia adelante

"✓ ◆2 ✓ ◆ #
z 1 z 1
+ a1 + a0 x(kT ) = b u(kT )
T T

Posterior al desarrollo de la expresión se llega a que

x(kT + 2T ) + (a1 T 2) x(kT + T ) + (1 a1 T + a0 T 2 ) x(kT ) = b T 2 u(kT )

Tarea: Utilizar los dos métodos restantes para obtener una versión discreta aproximada del sistema continuo
y comparar los resultados en simulación.

3. Transformada z
3.1. La Función de Transferencia Discreta
El objetivo de esta unidad es obtener la función de transferencia de un sistema lineal discreto invariante en
el tiempo por medio del análisis de la transformada z, donde z es una variable compleja.

9
 Considere el sistema discreto mostrado en la sig. Figura....
Si se aplica una secuencia de entrada u(tk ) = u(k T ) para k 2 Z+ [ 0, entonces el sistema generara una
salida y(kT ). Para que un sistema pueda ser descrito por ecuaciones de diferencia con coeficientes constantes,
debe ser lineal invariante en el tiempo.
Por ejemplo:
3y(k) + 2y(k 1) y(k 2) = 2u(k 1) 3u(k 2)
1
y(k) = [2u(k 1) 3u(k 2) 2y(k 1) + y(k 2)]
3
Si u(k) es una secuencia escalón unitario y las condiciones iniciales son y( 1) = 2, y( 2) = 1. entonces se
puede calcular y(k) de forma recursiva como: (Tabla...)
Una de las ventajas de este tipo de ecuaciones es que la solución se puede encontrar por sustitución directa,
sin embargo, la solución no esta en forma cerrada y es muy difícil extraer y/o analizar propiedades de la misma,
tales como estabilidad o diseño de un controlador.
La Transformada Z convierte una señal definida en el dominio del tiempo discreto en una representación
en el dominio de la frecuencia compleja, es decir, mapea de x(kT ) a X(z), donde z es una variable compleja.
La razón de ser de la transformada z es introducir una herramienta formal que permita extender la noción de
función de transferencia a los sistemas discretos y de esta forma tratar de manera unificada la parte continua y
la parte discreta en un sistema de control digital.

3.2. Definición de la Transformada z

Considere la siguiente señal muestreadora ideal
1
X
m(t) = (t kT )
k= 1
⇢
1 para t = kT
donde (t kT ) =
0 para t 6= kT
Sea f ⇤ (t) la salida de un muestreador ideal (Figura...), la cual está descrita por
1
X 1
X
f ⇤ (t) = f (t) (t kT ) = f (kT ) (t kT ) (8)
k=0 k=0

donde f (t) es la entrada al muestreador, comenzando el muestreo en el instante t = 0 [6].

Recuérdese que cuando hay un dispositivo de retención de orden cero después del muestreador de tipo pulso,
sólo importan los valores de la función f (t) en los instante de muestreo; la salida del retenedor es la misma,
sin interesar si la aproximación utilizada para el muestreador sea máximamente plana o ideal. Considerando
la transformada de Laplace del impulso desplazado un tiempo kT como e kT s = L{ (t kT )}, entonces la
transformada de Laplace de (8) se obtiene como
1
X
F ⇤ (s) = f (kT )e kT s

k=0

expresión que es la transformada de Laplace de la salida del muestreador ideal (o bien del retenedor de orden
cero).
Dado que la expresión F ⇤ (s) contiene el termino e T s , a diferencia de la mayoría de las FT de sistemas
continuos, la transformada F ⇤ (s) no es una función racional de s, por lo que se tendrán problemas al tomar la
transformada inversa de Laplace. De la necesidad de transformar la función irracional2 (dado que la exponencial
2 ex
P1 xn
= n=0
n!

10
es un numero irracional) F ⇤ (s) en una racional F (z), mediante una transformación de la variable compleja s en
otra z. Una selección para esta transformación es tomar

z = eT s .

O de forma inversa se puede obtener s = 1

T ln z. La relación entre s y z puede definirse como la transformación
z. entonces se tiene que

F (z) = F ⇤ (s)
s= T1 ln z
1
X
k
= f (kT )z
k=0

la cual cuando se escribe en forma cerrada, es una función de racional de z. Por consiguiente, puede definirse
F (z) como la transformada z de f (t), es decir, F (z) = Z[f (t)].
Si una señal x(tk ) que tiene valores discretos x0 , x1 , . . . , xk , . . ., entonces se define la transformada z de la
señal como

X(z) , Z [x(kT )]
1
X
k
= x(kT )z , t0 < |z| < R0 (9)
k= 1

donde y z es, en general, un número complejo de la forma z = A ej ! , donde A es el módulo de z, y ! es la

frecuencia angular en radianes por segundo (rad/s). Ademas, se asume que se pueden encontrar valores de r0
y R0 como límites sobre la magnitud de la variable compleja z tal que la serie (9) converge [4].
Observese que la expansión de (9) da como resultado
1 2 k
X(z) = x(0) + x(T )z + x(2T )z + · · · + x(kT )z + ··· (10)

La ecuación (10) implica que la transformada z de cualquier función en tiempo continuo x(t) se puede
escribir, mediante inspección, en la forma de serie. La z k en esta seria indica la posición en el tiempo en la que
se presenta la amplitud x(kT ). De manera contraria, si X(z) está dada en la forma de una seria como la que se
indicó, la trasformada z inversa se puede obtener por inspección como una secuencia de la función x(kT ) que
corresponde a los valores de x(t) en los valores de tiempo respectivos [9].

3.2.1. Región de Convergencia (ROC)

La región de convergencia, también conocida como ROC, define la región donde la transformada-z (TZ)
existe. La ROC es una región del plano complejo donde la TZ de una señal tiene una suma finita. La ROC para
una x(k) es definida como el rango de z para la cual la transformada-z converge. Ya que la transformada–z es
una serie de potencia, converge cuando x(k)z k es absolutamente sumable. Su definición formal es:
( 1
)
X
ROC = z : x(k)z k < 1 .
k= 1

Al utilizar el método de la transformada z para resolver problemas en tiempo discreto no es necesario

especificar los valores de z para los cuales X(z) converge.
Cuando en la ecuación (9) se permite que k pueda variar de 1 a 1, recibe el nombre de transformada
bilateral , mientras que para el caso en que k varia de 0 a 1, recibe el nombre de transformada unilateral
O bien, en [8].
Ts (T s)2 (T s)3
1 + + ···
e Ts
= 2 8 48
2 3
Ts (T s) (T s)
1+ + + + ···
2 8 48

11
 P1
X(z) = k=0 x(kT )z k , misma que para problemas reales es la que se utiliza al considerar que para valores de
k < 0, la señal x(kT ) = 0. La transformada unilateral, a menos que se indique lo contrario, será tratada en este
curso.
Por otro lado, la recuperación de la señal discreta x(kT ) a partir de de X(z), se realiza por medio de la
transformada Z inversa. Si la transformada z está dada como el cociente de los polinomios en z, entonces
la transformada z inversa se puede obtener por varios métodos diferentes, tales como el método de la división
directa, el método computacional, el método de expansión en fracciones parciales y el método de la integral de
inversión [9].
Ejemplo: Los datos x(kT ) son tomados como muestras de la señal x(t) = e at 1(t) con un periodo de
muestreo T , donde 1(t) es la función del escalón unitario, con valor cero para t < 0 y uno para t 0 (el
escalón unitario se agrega para definir el intervalo de tiempo para el cual esta definida la señal x(t)). Entonces
x(kT ) = e akT 1(kT ). Determine la transformada z de esta señal.
Solución. Considere la serie
1
X 1
1 + r + r2 + · · · = rk =
0
1 r
donde r es una constante real o compleja con magnitud menor de 1, esto es, para que la serie sea convergente
se requiere que |r| < 1, condición que es conocida también como el dominio de convergencia.
Aplicando (9), se tiene que
1
X 1
X
k akT k
x(kT )z = e z
k= 1 k=0
X1
aT 1 k
= e z
k=0
1
= aT z 1
1 e
z aT
= aT
, |z| > e , o bien, e aT
z 1
< 1.
z e
1 z
Ejemplo: Encuentre la transformada z del escalón con una amplitud dada A. Sol. U (z) = 1
= .
1 z z 1
zT
Tarea: Encuentre la transformada z de la función rampa f (t) = t 1(t). Sol. F (z) = 2 [2].
(z 1)
Tarea: Encuentre la transformada z de f (t) = sin wt. Sol. Considere las siguientes identidades: sin ✓ =
ei✓ e i✓ ei✓ + e i✓ z 1 sin wT
y cos ✓ = . Por lo tantoF (z) = .
2i 2 1 2z 1 cos wT + z 2
1 z 1 cos wT z 2 z cos wT
Tarea: Encuentre la transformada z de f (t) = cos wt. Sol. F (z) = =
1 2z 1 cos wT + z 2 z 2 2z cos wT + 1

3.2.2. Propiedades de la Transformada z

En esta sección se presentan las propiedades y teoremas mas relevantes de la transformada z que podrán ser
de utilidad al resolver problemas relacionados con la transformada z en aplicaciones de control.

Definición:
1
X
k
F (z) = f (kT ) z
k=0

Inversión [2]: I
1
f (kT ) = F (z) z k 1
dz
2⇡i

12
Linealidad [2]:
Z {↵ f (k) + g(k)} = ↵Z {f (k)} + Z {g(k)}
Prueba: A partir del lado derecho de la igualdad se tiene que
1
X
k
Z {↵ f + g} = (↵ f (k) + g(k)) z
k=0
1
X 1
X
k k
= ↵ f (k) z + g(k) z
k=0 k=0
= ↵Z {f (k)} + Z {g(k)}

k k
Ejemplo: Sea la función en tiempo discreto f (kT ) = 3 (2) 2 ( 1) , determine F (z). Sol. F (z) =
3 2
1 2z 1 1+z 1

Corrimiento en el Tiempo [2]:

n n
Z q f (k) = z F (z)
donde q es el factor de corrimiento.

Transformada de Convolución [2]:

( k ) ( k )
X X
Z {f ⇤ g} = Z f (n)g(k n) = Z f (k)g(n k) = (Z {f }) (Z {g}) = (Z {g}) (Z {f })
n=0 n=0

Teorema del Valor Inicial Si x(t) tiene la transformada z dada como X(z) y si el lı́mz!1 X(z) existe,
entonces el valor inicial x(0) de x(t) ó x(k) está dado por

x(0) = lı́m X(z) (11)

z!1

Para probar este teorema, observe que

1 2 k
X(z) = x(0) + x(T )z + x(2T )z + · · · + x(kT )z + ···

Al hacer que z ! 1en esta última ecuación, se obtiene la ecuación (11).

El teorema del valor inicial es conveniente para verificar la incidencia de posibles errores en el cálculo de la
transformada z. Debido a que x(0) se conoce, una verificación del valor inicial mediante lı́mz!1 X(z) puede
facilitar descubrir errores en X(z), en caso de existir.
Ejemplo: Determine el valor inicial x(0) si la transformada z de x(t) = 1 e t
está dada por (Nótese por
simple inspección que se espera obtener x(0) = 0)

1 e T z 1
X(z) =
(1 z 1 ) (1 e T z 1)

Mediante el uso del teorema del valor inicial, se encuentra que

1 e T z 1
X(z) = lı́m =0
z!1 (1 z 1 ) (1 e T z 1)

13
Teorema del Valor Final: suponga que x(k), donde x(k) = 0 para k = 0, tiene la transformada z dada
como X(z) y que todos lo polos de X(z) están dentro del círculo unitario, con la posible excepción de un solo
polo en z = 1 (Ésta el la condición para la estabilidad de X(z), o la condición para que x(k) (k = 0, 1, 2, ...)
permanezca finita.) Entonces el valor final de x(k), esto es, el valor de x(k) a medida que la k tiende a infinito,
puede darse mediante ⇥ ⇤
lı́m x(k) = lı́m 1 z 1 X(z)
k!1 z!1

Prueba: Considerando que

1
X
1 1 k
lı́m 1 z X(z) = lı́m 1 z x(kT )z
z!1 z!1
k=0
" 1 1
#
X X
k k
= lı́m x(kT )z x(kT T )z
z!1
k=0 k=0
" N N
!#
X X
k k
= lı́m lı́m x(kT )z x(kT T )z
z!1 N !1
k=0 k=0
" N N
!#
X X
k k
= lı́m lı́m x(kT )z x(kT T )z
N !1 z!1
k=0 k=0
" N N
#
X X
= lı́m x(kT ) x(kT T)
N !1
k=0 k=0

Desarrollando sumatorias
1
lı́m 1 z X(z) = lı́m [(x(0) + x(T ) + x(2T ) + · · · + x(N T T ) + x(N T )) (x( 1) + x(0) + x(T ) + x(2T ) + · · · +
z!1 N !1
= lı́m [x(N T )]
N !1
= x(1)

El teorema del valor final es muy útil para determinar el comportamiento de x(k) a medida que k ! 1 a
partir de su transformada z, esto es X(z).
Ejemplo: Determine el valor final de x(1) de
1 1
X(z) = 1 aT z 1
, a>0
1 z 1 e
mediante el uso del teorema del valor final. Note que la transformada X(z) es calculada a partir de la señalx(t) =
1 e a t , misma que evaluando x(1) = lı́mt!1 (1 e a t ) = 1.
Aplicando el teorema del valor final, se obtiene
⇥ ⇤
x(1) = lı́m 1 z 1 X(z)
z!1
 ✓ ◆
1 1
= lı́m 1 z 1
z!1 1 z 1 1 e aT z 1
✓ ◆
1 z 1
= lı́m 1
z!1 1 e aT z 1
= 1.

Ejemplo: Utilizando el Teorema de Valor Final y considerando la ROC correspondiente, determine x(1)
de
x(kT ) = ak 1(k).
Grafíque la evolución de x(kT ) para a = 0.8 y a = 0.8.

14
 Un resumen de las propiedades de la transformada z y algunas transformadas importantes las encontramos
en las siguiente tabla [9]

15
3.3. La Función de Transferencia de Pulso
En esta sección se introduce la transformada z como un mecanismo para la solución de ecuaciones en
diferencias. Las ecuaciones en diferencias se pueden solucionar fácilmente mediante el uso de una computadora,
siempre que se proporcionen los valores numéricos de todos los coeficientes y los parámetros iniciales. Sin
embargo, las expresiones en forma cerrada para x(k), usualmente no se pueden obtener a partir de la solución
por computadora. La utilidad del método de la transformada z es que permite obtener la expresión en forma
cerrada para x(k) [9].
Considere un sistema en tiempo discreto, lineal e invariante en el tiempo caracterizado por la siguiente
ecuación en diferencias:

a0 x(k + n) + a1 x(k + n 1) + · · · + an x(k) = b0 u(k + m) + b1 u(k + m 1) + · · · + bm u(k) (12)

16
donde u(k) y x(k) son la entrada y salida del sistema, respectivamente, en la k ésima iteración; los valores
a0 , a1 , ..., an . y b0 , b1 , ..., bm son constantes reales, donde n define el orden del sistema y donde ademas se
considera n m para que el sistema sea causal. Al describir la ecuación en diferencias en el plano z, se toma
la transformada z de cada uno de los términos en la ecuación, haciendo uso principalmente de los corrimientos
hacia atrás o hacia adelante, según corresponda. De esta forma, se tiene que la ecuación en diferencias puede
escribirse en el dominio de la variable compleja z como[9]

" n
# " #
X1 n
X2
n k n 1 k
a0 z X(z) x(k)z + a1 z X(z) x(k)z + · · · + an X(z) =
k=0 k=0
" m
# " #
X1 m
X2
b0 z m U (z) u(k)z k
+ b1 z m 1
U (z) u(k)z k
+ · · · + bm U (z)
k=0 k=0

Por facilidad de análisis, considere condiciones iniciales cero para la ecuación anterior, es decir, las sumatorias
representando las condiciones iniciales del sistema, son tomadas como cero, entonces se obtiene
a0 z n X(z) + a1 z n 1
X(z) + · · · + an X(z) = b0 z m U (z) + b1 z m 1
U (z) + · · · + bm U (z)
o bien

a0 z n + a1 z n 1
+ · · · + an X(z) = b0 z m + b1 z m 1
+ · · · + bm U (z)
que podemos escribir en forma de función de transferencia en términos de la variable compleja z como
X(z) b0 z m + b1 z m 1
+ · · · + bm
G(z) = = , m  n.
U (z) a0 z n + a1 z n 1 + · · · + an

Entonces podemos obtener una representación de la ecuación en diferencias (12) por medio de la multipli-
cación
X(z) = G(z) U (z)
misma que recibe el nombre de Función de Transferencia de Pulso o Función de Transferencia Dis-
creta del Sistema, que al considerar U (z) como un pulso, entonces U (z) = 1, entonces G(z) es igual a X(z)
ante una entrada tipo impulso.

Proponer un esquema donde se resuma e ilustro lo hasta ahora visto en esta unidad: el esquema
retro-alimentado de control, donde el controlador continuo puede ser representado por una ecua-
ción en diferencias, obtenido por una aproximación (Euler p.e.), y después obtener la función de
transferencia para tal ecuación en diferencias por medio de la transformada z.

3.3.1. Solución de la Ecuación de Estados Mediante la Transformada z

Ejemplo: Obtenga la función de transferencia pulso de la siguiente ecuación en diferencias empleando el

método de la transformada z:
x(k + 2) + 3 x(k + 1) + 2x(k) = 0, x(0) = 0, x(1) = 1.

Sol. Por el principio de linealidad, obtenemos la transformada z de cada termino como

Z{x(k + 2)} = z 2 X(z) z 2 x(0) z x(1)
2
= z X(z) z

Z{x(k + 1)} = zX(z) z x(0)

= z X(z)

17
 Z{x(k)} = X(z)
Representando la ecuación en diferencias por su transformada z, se obtiene
z z z z 1 1
X(z) = = = =
z2 + 3z + 2 (z + 1) (z + 2) z+1 z+2 1+z 1 1 + 2z 1

Adicionalmente, por simple inspección, se puede observar que a partir de cada término se puede concluir
que la señal en transformada X(z) corresponde en el tiempo discreto a la señal

x(k) = ( 1)k ( 2)k

haciendo uso de
1
X 1
1 + r + r2 + · · · = rk = .
0
1 r

Ejemplo: Determine la función de transferencia pulso de

y(kT ) + 0.5y(kT T ) + y(kT 2T ) = u(kT ) + u(kT T ).

1+z 1 z2 + z
Sol. G(z) = = .
1 + 0.5 z 1 + z 2 z 2 + 0.5 z + 1
Un análisis detallado de como obtener la función x(k) a partir de X(z) se da en la sección posterior corres-
pondiente a la transformada z inversa. Por el momento se concluye que a partir de la ecuación en diferencias,
ésta se puede resolver en forma cerrada a partir de la transformada z y posteriormente aplicando la transformada
z inversa.
De manera general, se puede emplear la transformada z para resolver la Ecuación de Estado [2]

x(k + 1) = Ax(k) + Bu(k)

y(k) = Cx(k) + Du(k)

Pre-multiplicando ambos lados de la ecuación de estados y tomando la transformada z termino a termino se

obtiene !
X1 X1 1
X X1
k k
z x(k + 1) = z z x(k) x(0) = Az k x(k) + Bz k u(k).
k=0 k=0 k=0 k=0

Así

z (X(z) x(0)) = A X(z) + B U (z)

1
X(z) = (zI A) (z x(0) + B U (z))

y ⇣ ⌘
1 1
Y (z) = C (zI A) zx(0) + C (zI A) B + D U (z).

A partir de la ecuación anterior, se puede introducir la función de transferencia pulso como

1
G(z) = C (zI A) B + D.

La secuencia en el tiempo de y(k) puede obtenerse usando la transformada z inversa.

3.4. Polos y Ceros de la Función de Transferencia Pulso

En aplicaciones de ingeniería del método de la transformada z, G(z) puede tener la forma

b0 z m + b1 z m 1
+ · · · + bm
G(z) = , mn
a0 z n + a1 z n 1 + · · · + an

18
o bien, expresando la ecuación anterior en la forma de ceros y polos, considerando como ceros aquellas raíces
de G(z) tal que G(z) = 0 y como polo aquellas raíces que hacen a G(z) = 1, se puede escribir entonces

X(z) K (z z1 ) (z z2 ) · · · (z zm )
G(z) = =
U (z) (z p1 ) (z p2 ) · · · (z pn )
donde los pi , con i = 1, 2, . . . , n son los polos de G(z), mientras que los zj , con j = 1, 2, . . . , m son los ceros de
G(z).
De la función de transferencia de pulso, la ubicación de los polos y los ceros de G(z) determinan las caracte-
rísticas de x(k), la la secuencia de valores o números. Como en el caso del análisis de sistemas de control lineales
en tiempo continuo en el plano s, también se utiliza una representación gráfica de las localizaciones de los polos
y ceros de G(z) en el plano z.
En ingeniería de control y en procesamiento de señales, una función de transferencia G(z) a menudo se
expresa como un cociente de polinomios en z 1 , como sigue
(n m) (n m+1) n
b0 z + b1 z + · · · + bm z
G(z) = 1
a0 + a1 z + · · · + an z n

donde z 1 se interpreta como el operador retraso unitario. Entonces, ha de utilizarse la presentación tanto en
potencias positivas, o bien negativas, dependiendo de las circunstancias y conveniencia. Por ejemplo, la función
de transferencia de pulso
z 2 + 0.5 z z (z + 0.5)
G(z) = 2 =
z + 3z + 2 (z + 1) (z + 2)
cuya ubicación de los polos esta en z = 1yz= 2, mientras que los ceros en z = 0 y z = 0.5, o bien

1 + 0.5 z 1 1 + 0.5 z 1
G(z) = =
1 + 3z 1 + 2z 2 (1 + z 1 ) (1 + 2 z 1)

cuya ubicación de los polos esta en z = 1 y z = 2, mientras que solo aparece el cero z = 0.5 a simple vista
y el cero ubicado en el origen desaparece en esta representación de manera explícita, por lo que es recomendable
para fines de ubicación de los polos y ceros, la representación en potencias positivas.

3.4.1. Raíces Dominantes de la Ecuación Característica

En general, algunas de las raíces, dada la ubicación de éstas en el plano z, afectan más la respuesta del
sistema que otras. Para fines de análisis y diseño, es importante aislar las raíces con un efecto dominante sobre
la respuesta transitoria y denominarlas raíces dominantes.
En el plano s, las raíces más próximas al eje j! en el semiplano izquierdo son las raíces dominantes [6],
ya que corresponden a respuestas en el tiempo que decaen con lentitud. Algunos autores distinguen los polos
dominantes de los no dominantes si se considera una separación entre ellos de al menos 10 veces, es decir, D = 10
en la figura siguiente:

19
 Hablando de un sistema que tiene polos con parte real solamente s = , la respuesta en el tiempo seria
de la forma e t mientras que para un polo complejo s = ± j! la respuesta en el tiempo sería de la forma
e t cos (!t) [3]. De ahí que para considerar si son los polos dominantes o no, solo importa la parte real con
respecto al eje j!, ya que esta es la que aparece en la exponencial decreciente para el caso estable. En [8] se
habla también de los polos dominantes en lazo cerrado.
De hecho, en [6] se puede obtener una reducción de un sistema de cuarto orden a uno de orden 2 considerando
solamente los polos dominantes, así una
K
G(s) =
(s + p1 ) (s + p2 ) (s2 + 2⇣!n s + !n2 )

se puede aproximar por

K
G(s) ⇠
=
p1 p2 (s2 + 2⇣!n s + !n2 )
en el caso de que p1 y p2 sean al menos seis veces mas grandes que !n .

3.5. Expresión alternativa para obtener F (z)

En el análisis de sistemas lineales es frecuente que la función de transferencia F (s) ya esté dada, y lo que
tenga que determinarse sea la transformada z, F (z), directamente a partir de F (s), sin tener que plantear
primero f (t) para después F (z) [6]. Para ello considere los siguientes casos.

4. Modelado desde el Punto de Vista del Proceso

4.1. Introducción
Hay básicamente dos formas de ver el esquema de control clásico de control digital en lazo cerrado como el
de la siguiente figura:

20
 1. Desde el punto de vista del proceso continuo: En este esquema el proceso continuo considera que el
algoritmo de control también es continuo dado que recibe una señal de control continua u(t) que proviene
del convertidor D/A y ademas de que su salida y(t) es de naturaleza continua, por lo que pasa desapercibido
que el control es realizado por un dispositivo digital, utilizando por ejemplo ecuaciones en diferencias, etc.
Con la idea de tener un esquema de control homogéneo (solo componentes continuos), surge la idea de
ver por las implicaciones que tendría el modelar el sistema de control discreto por una representación
continua, sin embargo esta idea tiene un inconveniente importante, que es el hecho de que entre los
instantes de muestreo no se tiene información, y el esquema resultante tendrá errores significativos debido
a la conversión de digital a continuo. De este hecho, se plantea en su lugar, convertir los elementos
continuos a su versión discreta, en el entendido que la información en los instantes de muestreo entre la
parte continua y la discreta coinciden.
2. Desde el punto de vista del elemento discreto (computadora): Debido a que el elemento discreto, corres-
pondiente al algoritmo de control, ve al sistema continuo como si fuese discreto (dado que este recibe una
señal digital y(k) y su salida es digital u(k) por naturaleza. Resulta conveniente por lo tanto modelar el
sistema continuo por uno discreto.

Por medio de la transformada z se puede representar lo que ocurre dentro de la computadora en términos
de su función de transferencia de pulso, de manera que el esquema general de control digital o el algoritmo
de control, se puede re-plantear en términos de funciones de transferencia. Nótese que aunque se cuente con
funciones de transferencia tanto para la parte discreta como para la continua, ambas se encuentran en distintos
dominios. Con el fin de contar con un modelo homogéneo, el convertidor de analógico a digital se modelado por
un muestreador simplemente, mientras que el convertidor de digital a analógico se modela por un muestreador
seguido de un reconstructor, siendo el ZOH (Retenedor de orden cero) el mas utilizado y por ende el que en
este curso se estudiará.
Para la obtención del modelo homogéneo, se procede a obtener la función de transferencia z equivalente del
retenedor de orden cero en cascada con G(s), es decir, se parte de u(k) en la figura siguiente y se obtiene el
modelo matemático hasta la salida ficticia y(k), que está al final de la salida del proceso continuo (véase la
figura siguiente).

La salida ficticia y(k) se introduce por dos razones principales:

1. Obtener una función de transferencia G(z) que relaciones la salida continua y(t) con la entrada u(k), esto
para contar con una salida discreta y(k) y poder obtener G(z) con entrada u(k).
2. Representar la perdida de información entre los instantes de muestreo.

21
4.2. Obtención de G(z) con retenedor de orden cero a partir de G(s)
El primer paso para realizar el diseño y análisis del sistema de control digital, el cual contiene elementos
continuos, es encontrar la función de transferencia discreta de la porción continua. Para un sistema como el que
se muestra en la figura anterior, se desea encontrar la función de transferencia entre u(k) y y(k). Ademas, existe
una equivalencia exacta discreta para este análisis debido a que la retención de orden cero describe precisamente
lo que sucede entre las muestras, mientras que la salida y(k) solamente depende de la entrada en los tiempos
de muestreo u(k) [3].
El efecto del retenedor de orden cero puede verse como la acción de dos escalones

b0 = u(t) u(t T ), t = kT, k = {0, 1, 2, · · · }

⌧s ⌧s
1 e 1 e
b0 (s) = =
s s s
(La ultima utilizando la propiedad de la transformada de Laplace de un escalón con retardo)
Para una planta descrita por una G(s) y precedida por un retenedor de orden cero, la función de
trasferencia discreta es ⇢
G(s)
G(z) = 1 z 1 Z (13)
s
considerando el cambio de variable z = e⌧ s . La formula tiene un termino G(s)/s debido el control aparece
como una entrada escalón durante cada periodo de muestreo. Se tiene el término 1 z 1 debido a que un
escalón de duración de una muestra se puede conceptuar como un escalón de duración infinita seguido de un
escalón negativo (o bien positivo según el caso) con retraso de un ciclo. La formula de la ecuación (13) permite
reemplazar el sistema mixto (continuo y discreto) que se muestra en la Fig. 8.12a por el sistema equivalente
discreto puro que se muestra en la Fig. 8.12(b).

El análisis y diseño de sistemas discretos es muy similar al análisis y diseño de sistemas continuos; de hecho,
se aplican las mismas reglas. La función de trasferencia de lazo cerrado de la Fig. 8.12b se obtiene empleando

22
las mismas reglas de la reducción del diagrama de bloques, esto es

Y (z) D(z)G(z)
= (14)
R(z) 1 + D(z)G(z)

Como se requiere encontrar el comportamiento característico del sistema de lazo cerrado, se desea encontrar los
factores del denominador de la ecuación (14); esto es, encontrar las raíces reales de la ecuación característica
([3].)
1 + D(z)G(z) = 0.

4.2.1. Reducción de Bloques

Cuando se tienen sistemas de control digital en donde intervienen tanto partes continuas, discretas y bloques
de discretización, es conveniente realizar reduccion para facilidad de análisis y diseño. De esta forma considere

Y (z)
donde G(z) = Z {G(s)} y H(z) = Z {H(s)}. Por lo tanto = G(z)H(z).
U (z)
Ahora considere el caso sin muestreador intermedio

Y (z)
donde HG(z) = Z {G(s)H(s)}. Por lo tanto = HG(z).
U (z)

23
4.3. Discretización de funciones de transferencia en lazo abierto usando tablas de
transformadas
El procedimiento para transformar una expresión del dominio de Laplace al dominio z, se puede realizar al
transformar de del dominio s al tiempo y posteriormente del tiempo al dominio z, o bien, a partir de tablas de
transformadas. Para realizar esto tenemos que partir de expresiones de G(s) en forma de función racional de s
y aplicar el procedimiento de expansión en fracciones parciales, las cuales son fáciles de encontrar en las tablas,
como se ilustra en los siguientes ejemplos.

4.3.1. Ejemplo: Sistemas de Primer Orden

Determine la función de transferencia G(z) a partir de la función de transferencia G(s) considerando que
esta última es precedida por un retenedor de orden cero, con
b
G(s) =
s+a

Sol. La función de transferencia se determina como

⇢
b
G(z) = 1 z 1 Z
s (s + a)
⇢1 1
1
= b 1 z Z a a
s (s + a)
⇢ ⇢ ⇢
1 1 1 1 1
= b 1 z Z Z
a s a s+a
⇢  
1 1 1 1
= b 1 z 1
a 1 z 1 a 1 e aTz 1

b 1 z 1
= 1
a 1 e aTz 1
" #
b 1 e aTz 1 1 z 1
=
a 1 e aTz 1
b 1 e aT z 1
=
a ae aTz 1

b 1 e aT
=
a z e aT
Y (z)
=
U (z)

Al representar la función de transferencia anterior como una ecuación en diferencias se obtiene que

aT b aT
y [(k + 1)T ] = e y [kT ] + 1 e u [kT ]
a
mismo que es el resultado de la discretización exacta de un sistema de primer orden, visto con anterioridad.
Tarea: Obtenga la representación discreta de un sistema de segundo orden que es precedido por un retenedor
de orden cero, donde
1
G(s) =
s (s + 3)
Tarea: Usando tablas, obtenga la representación discreta la siguiente función de transferencia que es prece-
dida por un retenedor de orden cero, donde
1
G(s) = 2
s
El código de Matlab para su obtención es el siguiente considerando T = 1 s:

24
T = 1;
numC=1;
denC=[1 0 0];
sysC=tf(numC,denC);
sysD=c2d(sysC,T);
[numD,denD,T] = tfdata(sysD)

cuyo resultado da
numD = [0 0.5 0.5] y denD = [1 2 1]
0.5z + 0.5 z+1
que corresponde a G(z) = 2 = 0.5 2.
z 2z + 1 (z 1)
Ejemplo: Obtenga la representación discreta del siguiente sistema que es precedido por un retenedor de
orden cero
1
G(s) = 2
s +s+1
Tarea: Obtenga la representación discreta de un sistema de segundo orden que es precedido por un retenedor
de orden cero, donde
s+3
G(s) = 4
s + 4s + 6s2 + 5s + 2
3

4.3.2. Modelado de un Motor de CD

Considere el modelado de un motor de corriente directa.

Este motor provee de un movimiento rotatorio o bien, conectado a bandas transportadoras, cables, etc., tal
que se produzca un movimiento traslacional.
El torque que puede producirse por un motor viene directamente proporcional a la corriente de armadura
(corriente del rotor), el cual puede ser descrito por:

⌧ = Kt i.

La fuerza electromotriz (emf), e, esta relacionada a la velocidad rotacional por la siguiente expresión:

˙
e = Ke ✓.

Para las expresiones anteriores Kt = Ke , misma que por generalidad llamaremos simplemente K.
De las leyes de Newton en combinación con las de Kirchhoff se obtiene que:

J ✓¨ + b ✓˙ = ⌧ = K i

di di
L + R i + e = L + R i + K ✓˙ = V.
dt dt

25
donde el termino J ✓¨ es el equivalente a masa por aceleración en la segunda ley de Newton y b ✓˙ representa la
fuerza de fricción entre el rotor y el estator, etc.
La
Tarea: Determinar la función de transferencia en z que represente la relación entrada-salida, considerando
el voltaje aplicado como la variable de entrada y la velocidad angular como la salida a controlar. Utilice un
periodo de muestreo de 0.01 s ó 0.001 s. Para el mismo sistema del motor, determine su representación en
espacio de estados.
Simular ambas respuestas ante un voltaje de entrada unitario.
˙
✓(s) K 2
= = 2
V (s) (sL + R) (sJ + b) + K 2 s + 12 s + 20.02

para los parámetros del modelo dados como J = 0.01, b = 0.1, K = 0.01, R = 1, L = 0.5. La obtención de la
función de transferencia en tiempo discreto es dada a partir del siguiente código de Matlab que resulta en (con
un periodo de muestreo de 0.12 segundos)

˙
✓(z) 0.0092 z + 0.0057
= 2
V (z) z 1.0877z + 0.2369

El código de Matlab para su obtención es el siguiente:

R=1;
L=0.5;
Kt=0.01;
J=0.01;
b=0.1;
num = Kt;
den = [(J*L) (J*R)+(L*b) (R*b)+(Kt^2)];
Ts = 0.12;
[numz,denz] = c2dm(num,den,Ts,’zoh’)
[x1] = dstep(numz,denz,101);
t=0:0.12:12; stairs(t,x1)
xlabel(’Time (seconds)’)
ylabel(’Velocity (rad/s)’)

4.4. Escribir la sección 4.3 del libro [4], pp. 59.

Muy útil el ejemplo 4.5 del mismo libro en la pp. 60, sobre la obtención del modelo discreto a partir del
continuo considerando que es precedido por el ZOH.

4.5. Transformada Z inversa

Se han establecido dos formalismos para la descripción de los sistemas discretos: la ecuación en diferencias y
la función de transferencia. En esta sección se aborda el análisis del comportamiento o respuesta de los sistema
discretos a partir de determinadas condiciones iniciales y/o entradas de prueba. Existen diferentes métodos para
dar solución a la ecuación en diferencias o a la función de transferencia en z, tales como

1. Soluciones iterativas;
2. Transformada Z inversa;
3. Métodos analíticos.

26
4.6. Solución a una ecuación en diferencias por iteración numérica
Este es un procedimiento sencillo que consiste en determinar la variable de salida de un sistema a partir de
las condiciones iniciales sustituidas en la ecuación en diferencia, pero tiene la desventaja que es difícil determinar
el comportamiento del sistema de forma cerrada.

4.6.1. Ejemplo: sistema de primer orden

Considere la ecuación en diferencias correspondiente a un sistema de primer orden ante una entrada escalón
y con condiciones iniciales cero dada como
y [(k + 1)T ] = ↵y(kT ) + u(kT ), y(0) = 0

Recordar el programa de excel....

Analizando las operaciones realizadas por el procedimiento iterativo obtendríamos:
k =0 y(T ) = ↵y(0) + u(0) = ↵ · 0 + · 1 =
k =1 y(2T ) = ↵y(T ) + u(T ) = ↵ · + · 1 = (↵ + 1) ·
k =2 y(3T ) = ↵y(2T ) + u(2T ) = ↵ · (↵ + 1) · + · 1 = ↵2 + ↵ + 1 ·
k =3 y(4T ) = ↵y(3T ) + u(3T ) = ↵ · ↵2 + ↵ + 1 · + · 1 = ↵3 + ↵2 + ↵ + 1 ·
de la cual se puede deducir que para k = n 1:
y(nT ) = ↵y(nT T ) + u(nT T ) = ↵n 1
+ ··· + ↵ + 1 ·

Recordando la sumatoria geométrica

1 ↵n
↵n 1
+ ··· + ↵ + 1 =
1 ↵
Por lo que
1 ↵n
y(nT ) =
1 ↵
Considerando T = 1, = 1 y analizando los siguientes casos para valores de ↵ = 0.9, 0.9, 1.1111, 1.1111.
Ver gráficas...
Como puede observarse va a depender del valor que adopte ↵ para obtener un comportamiento convergente
a un valor (comportamiento estable) o bien divergente hacia el infinito (comportamiento inestable).
La obtención de la solución para y(kT ) obtenida de forma iterativa, dependerá por mucho del orden del
sistema que se este analizando, dando lugar a que resulte bastante complejo la obtención de una expresión
analíticas como la anterior para un sistema de orden mayor.

4.7. Solución de la Ecuación Lineal de Diferencias Mediante Transformada Z In-

versa
Se puede determinar la respuesta de un sistema con salida y(kT ), ante una entrada u(kT ) y bajo condiciones
iniciales cero, a partir de su función de transferencia
Y (z)
G(z) = .
U (z)

La salida del sistema y(kT ) se puede determinar a partir de la función de transferencia como
Y (z) = G(z)U (z)
y finalmente por medio de la transformada Z inversa definida como
1
y(kT ) = Z {Y (z)}
misma que puede determinarse por tablas o por manipulaciones algebraicas a fin de obtener y(k) a partir
expresiones sencillas que tengan transformada Z inversa directa desde Y (z).

27
4.8. Cálculo Analítico de la Transformada Z Inversa
Al proceso de obtener la señal en el tiempo discreto y(kT ) a partir de su expresión transformada Y (z) se le
denomina transformada Z inversa. Esta herramienta juega el mismo papel que la transformada de Laplace en
los sistemas de control en tiempo continuo. Se hace notar que la transformada z de inversa da como resultado
una única y(k). Esto significa que la transformada inversa dará como resultado una secuencia de tiempo que
especifica los valores de y(t) solamente en los instantes de tiempo discreto.
Ademas de utilizar tablas para la determinación de la transformada z inversa, podemos enunciar cuatro
métodos que dan solución a la transformada z inversa:

1. Método de expansión en fracciones parciales (resultado en forma cerrada);

2. Método de la división directa (resultado en forma de series);
3. Método computacional (resultado en forma de series);
4. Método de la integral de inversión (resultado en forma cerrada).

4.8.1. Método de Fracciones Parciales

Aun con una tabla extensa de formulas de la transformada z inversa, se encontrará dificultad para determinar
la señal y(k) si el sistema bajo estudio es una función complicada. En estos casos, como se lleva a cabo en los
sistemas continuos, se utiliza en método de la expansión en fracciones parciales con la finalidad de llegar a
expresiones que se puedan reconocer fácilmente a partir de tablas de la transformada z. Una de las grandes
ventajas de este método es que proporciona una solución en forma cerrada.
Y (z)
El procedimiento consiste en determinar la expansión en fracciones parciales de , donde la división por
z
z se debe a que es usual encontrar este termino en las transformadas z en el numerador, por ejemplo:
z
Función escalón
z 1
T z2
Función rampa 2
(z 1)

Ejemplo: Obtenga la función de transferencia pulso de la siguiente ecuación en diferencias empleando el método
de la transformada z:

x(k + 2) + 3 x(k + 1) + 2x(k) = 0, x(0) = 0, x(1) = 1.

Sol. Por el principio de linealidad, obtenemos la transformada z de cada termino como

Z{x(k + 2)} = z 2 X(z) z 2 x(0) z x(1)

2
= z X(z) z

Z{x(k + 1)} = zX(z) z x(0)

= z X(z)

Z{x(k)} = X(z).
Representando la ecuación en diferencias por su transformada z, se obtiene
z z z z 1 1
X(z) = = = = .
z2 + 3 z + 2 (z + 1) (z + 2) z+1 z+2 1+z 1 1 + 2z 1

28
 Adicionalmente, por simple inspección, se puede observar que a partir de cada término se puede concluir
que la señal en transformada X(z) corresponde en el tiempo discreto a

x(k) = ( 1)k ( 2)k

haciendo uso de
1
X
2 1
1 + r + r + ··· = rk = .
0
1 r

Hacer un comentario de que si deseáramos tener el valor de y(kT ) con k = 1000 tendríamos que evaluar 50
veces la ecuación en diferencias, mientras que con la solución cerrada basta con calcular y(k) a partir del valor
de k correspondiente..
Tarea: evaluar en excel o Matlab la solución de la ecuación en diferencias y la solución por fracciones
parciales.
Ejemplo: Determine la transformada z inversa de
aT
1 e z
Y (z) = aT )
.
(z 1) (z e

Tarea: determine la transformada inversa de

2z 3 + z
Y (x) = 2 .
(z 2) (z 1)

Ejemplo: Determinar la respuesta de un sistema de primer orden por medio de la transformada z inversa.
Sea el sistema discreto
y [(k + 1) T ] = ↵y(kT ) + u(kT )
donde ↵ y son constantes, la condición inicial y(0) = 0 y una entrada tipo escalón unitario u(kT ) = 1 para
k 0 y u(kT ) = 0 para k < 0.
Ejemplo: Polos Repetidos
Determine la respuesta en el tiempo discreto de la función de transferencia dada como
1
Y (z) = .
(1 + z 1 ) (1 z 1 )2

Expresando la anterior como una FT con potencias positivas

1
Y (z) =
1 1)
1 2 1 )2
z (1 + z z (1 z
z z2
z3
= 2.
(z + 1) (z 1)

Primeramente reescribimos la expresión anterior como:

Y (z) z2
= 2.
z (z + 1) (z 1)

Note que hay polos repetidos. Para este caso se sigue la siguiente forma:
A1 A2 An
+ 2 + ··· + n
z p (z p) (z p)
con 
1 dn i
n Y (z)
Ai = (z p) , i = 1, 2, . . . , n.
(n i)! dz n i z z=p

29
Se puede escribir entonces
Y (z) A1 A2 B
= + + .
z z p (z p)2 z+1
Los coeficientes son entonces (usando división del cociente3 ):

d 2 Y (z) 3
i = 1, A1 = (z 1) =
dz z z=1 4

2 Y (z) 1
i = 2, A2 = (z 1) =
z z=1 2
" #
z2 1
B= 2 = .
(z 1) z= 1 4

Entonces se tiene que

3 z 1 z 1 z
Y (z) = + + .
4 z 1 2 (z 1)2 4z+1

De las tabla se tiene que

3 1 1 k
y(k) = + k + ( 1) .
4 2 4

Tarea: Determine la respuesta en tiempo discreto de la siguiente función de transferencia

1+z 1
Y (z) =
1 z 1 + 0.5z 2

4.8.2. Método de la División Directa

En este método la transformada z inversa se obtiene mediante la expansión de Y (z) en una serie infinita
de potencias de z 1 y es utilizado cuando no es sencillo obtener una solución en forma cerrada o solo se esta
interesado en conocer los primeros términos de la señal y(k). El método proviene del hecho de que si X(z) está

3 f (x) u u0 v uv 0
= , entonces f 0 (x) = .
v v2

30
expandida en una serie de potencias de z 1
como
1
X
k
Y (z) = y(kT )z
k=0
1 2 k
= y(0) + y(T )z + y(2T )z + · · · + y(kT )z + ···

entonces y(kT ) es el coeficiente del termino z k

. Por los tanto, los valores de y(kT ) para k = 0, 1, 2... se pueden
determinar por inspección.
Ejemplo: Encuentre y(k) para k = 1, 2, 3, 4, cuando Y (z) esta dada por
10z + 5
Y (z) =
(z 1) (z 0.2)

Sol. Expresar Y (z) como polinomios en z k

, por lo tanto

10z 1 + 5z 2
Y (z) =
1 1.2z 1 + 0.2z 2

Al dividir el numerador entre el denominador, se obtiene que

De este modo
1 2 3
Y (z) = 10z + 17z + 18.68z 4 + · · ·
+ 18.4z
P1
Al compara esta expansión de Y (z) en una serie infinita con Y (z) = k=0 y(k)z k
, se obtiene que

y(0) = 0
y(1) = 10
y(2) = 17
y(3) = 18.4
y(4) = 18.68
..
.

En general no es posible, salvo en casos especiales, determinar una solución cerrada.

Tarea: Usando el método de la división directa, determine la transformada z de

z2
Y (z) =
z2 1.5z + 0.5
para k = 0, 1, 2, 3, 4
⇢
3 7 15
Sol. y(k) = 1, , , , ... .
2 4 8

31
4.8.3. Método de la Inversión Directa

Es una técnica utilizada para la obtención de la transformada z inversa mediante la integral de inversión
dada como I
1 1
y(k) = Z {Y (z)} = Y (z)z k 1 dz
2⇡j C
donde C es un circulo con centro en el origen del plano z tal que todos los polos de Y (z)z k 1
están dentro de
él [9].
La ecuación que resulta en la transformada z inversa en términos de los residuos se puede obtener si se
utiliza la teoría de la variable compleja. Ésta se puede obtener como
y(k) = K 1 + K 2 + · · · + Km
m
X ⇥ ⇤
= residuo de Y (z)z k 1
en el polo z = zi (15)
i=1

donde K1 , K2 , ... , Km denotan los residuos de Y (z)z k 1 en los polos z1 , z2 , ..., zm , respectivamente. Al evaluar
los residuos, observe que si el denominador Y (z)z k 1 contiene un polo simple en z = zi entonces el residuo K
correspondiente está dado por ⇥ ⇤
Ki = lı́m (z zi ) Y (z)z k 1 .
z!zi

Si Y (z)z k 1
contiene un polo múltiple zj de grado q, entonces el residuo K está dado por [9]
1 dq 1 ⇥ q ⇤
K= lı́m (z zi ) Y (z)z k 1
.
(q 1)! z!zi dz q 1

Un aspecto importante de este método con respecto a fracciones parciales, es que cuando se tiene un término
en el denominador con polos repetidos, solo se debe calcular un ÚNICO valor de K correspondiente al término
A1
de las raíces repetidas, mientras que en fracciones parciales hay que descomponer, por ejemplo, como 1 +
s
A2 A3 An
+ 3 + ··· + n .
s2 s s
Debe observarse que el método de la integral de inversión, cuando se evalúa por residuos, es una técnica muy
sencilla para obtener la transformada z inversa, siempre que Y (z)z k 1 , no tenga polos en el origen, dado que
esto se podría tornar tedioso. Para tales casos se recomienda el método de expansión en fracciones parciales.
Ejemplo: Obtenga y(k) por el método de la integral de inversión cuando Y (z) está dada por
z 1 e aT
Y (z) = aT )
.
(z 1) (z e
Multiplicando ambos lados de la expresión anterior por z k 1
se tiene que
1 e aT z k
Y (z)z k 1
= .
(z 1) (z e aT )
Para k = 0, 1, 2, ..., la expresión de Y (z)z k 1 tiene dos polos simples en z1 = 1 y z2 = e aT
. Por lo tanto a
partir de (15), se tiene
2
" #
X 1 e aT z k
y(k) = residuo de en el polo z = zi
i=1
(z 1) (z e aT )
= K 1 + K2
donde
K1 = [residuo en el polo simple z = 1]
" #
1 e aT z k
= lı́m (z 1) =1
z!1 (z 1) (z e aT )

32
 ⇥ ⇤
K2 = residuo en el polo simple z = e aT
" #
aT k
1 e z
= lı́m z e aT = e akT
z!e aT (z 1) (z e aT )

Por lo tanto
akT
y(k) = K1 + K2 = 1 e , k = 0, 1, 2, ...
a
que correspondería a la respuesta de un sistema en tiempo continuo de la forma G(s) = ante una
s+a
a
entrada tipo escalón unitario (Y (s) = ). Note que en este procedimiento no hay que aplicar tablas de
s (s + a)
transformadas inversas.
Ejemplo: Polos Repetidos
Mediante el método de la inversión directa, determine la respuesta en el tiempo discreto de la función de
transferencia dada como
1
Y (z) = 2.
(1 + z 1 ) (1 z 1 )

Sol.: Después de manipulaciones algebraicas se tiene que

z k+2
Y (z)z k 1
= 2.
(z + 1) (z 1)

Determinando el primer residuo para el polo simple se tiene

" #
z k+2
K1 = lı́m (z + 1) 2
z! 1 (z + 1) (z 1)
" #
z k+2
= lı́m 2
z! 1 (z 1)
k+2
( 1)
= 2
( 1 1)
k 2
( 1) ( 1)
= 2
( 1 1)
k
( 1)
= .
4
El siguiente residuo se calcula como

1 d2 1 h 2
i
K2 = lı́m (z 1) Y (z)z k 1
(2 1)! z!1 dz 2 1

d z k+2
= lı́m
z!1 dz z + 1
" #
(z + 1) (2 + k) z k+1 z k+2
= lı́m 2
z!1 (z + 1)
3 k
= + .
4 2
Finalmente

y(k) = K 1 + K2
k
3 k ( 1)
= + + .
4 2 4

33
4.9. Transformada Z Modificada
El comportamiento entre puntos de muestreo pueden ser investigados usando la transformada Z modificada.
En términos generales, la transformada Z modificada es la transformada ordinaria añadiendo un tiempo de
adelanto mT , el cual es una fracción del periodo de muestreo. Al modificar el tiempo de retardo mT es posible
recuperar toda la información que se quiera sobre la señal entre los instantes de muestreo[6]. La transformada
Z modificada se define como [2]
1
X
k
F (z, m) = z f (kT T + mT ), 0  m  1.
k=0

La siguiente figura ilustra el proceso de la transformada z mediante los siguientes tres pasos:

1. La función f (t) se desplaza mT unidades a la izquierda (adelanto en el tiempo), donde 0 < m < 1. Esto
resulta en f (t + mT ).
2. La función f (t + mT ) se muestrea con un muestreador ideal; este proceso se inicia en t = 0.
3. La secuencia del muestreador se recorre a la derecha (retraso en el tiempo) un instante de muestreo T .

Ejemplo: Determine la transformada z modificada de la función f (t) = e at

1(t) [6].
Sol. La transformada z modificada resulta en
1
X
a(kT T +mT ) k
F (z, m) = e z
k=0
1
X
1 amT a(kT ) k
= z e e z
k=0
amT
1 amT z e
= z e aT
= aT
.
z e z e

5. Análisis de Sistemas de Control en Tiempo Discreto

5.1. Introducción
En esta unidad se lleva a cabo un análisis sobre las propiedades y conceptos básicos involucrados en los
sistemas de control, tales como estabilidad, controlabilidad, observabilidad, etc.
Uno de los requisitos mas importantes en el desempeño de los sistemas de control es la estabilidad, tanto
para sistemas continuos como para discretos. Para tal fin considere el sistema lineal discreto

x(k + 1) = Ax(k) + Bu(k) (16)

y(k) = Cx(k) + Du(k)

Definición de Estabilidad: El punto de equilibrio xe = 0 de (16) es:

Estable, si para cada ✏ > 0 existe = (✏) tal que

kx(0)k < ) kx(t)k < ✏, 8t 0.

Inestable, en caso contrario

Asintóticamente estable (AE), si es estable y puede elegirse tal que

kx(0)k < ) lı́mt!1 x(t) = 0.

34
Estabilidad Entrada Acotada-Salida Acotada: El sistema descrito por (16) es estable en el sentido entrada
acotada-salida acotada (BIBO, bounded input bounded output), si para cualquier entrada acotada u(k), la salida
y(k) está acotada.
Teorema (Estabilidad y Raíces de la Ecuación Característica): Para que el sistema (16) sea asintóti-
camente estable, se requiere que las raíces de la ecuación característica se encuentren dentro del círculo unitario
en el plano complejo z.
Prueba: [6].
Por razones practicas, cuando la ecuación característica tenga al menos una raíz sobre el circulo unitario, se
dirá que el sistema es marginalmente estable o marginalmente inestable.

5.2. Análisis de Estabilidad de Sistemas Discretos

Para los sistemas discretos, la estabilidad se determina por medio del análisis de las raíces del polinomio
característico y que éstas estén dentro del círculo unitario en el plano z. El criterio de Routh-Hurwitz, el criterio
de Nyquist, el diagrama del lugar geométrico de las raíces y el diagrama de Bode, originalmente establecidos
para sistemas continuos, se pueden aplicar al estudio de la estabilidad de un sistema lineal en tiempo discreto.

5.2.1. Criterio de Routh-Hurwitz

Una de las herramientas útiles para probar estabilidad es el criterio de Routh-Hurwitz, cuyo principio se
basa en encontrar una transformación que modifique el plano complejo z en otro plano complejo. Una de las
transformaciones más simples es la transformada r (transformada bilineal) dada como:
1+r
z= .
1 r

Una vez que la ecuación característica en z se transforma al dominio r, se puede aplicar el criterio de
Routh-Hurwitz, como usualmente se hace para sistemas en tiempo continuo a la nueva ecuación en términos de
r.
De forma general, dado el polinomio característico:

P (z) = a0 z n + a1 z n 1
+ · · · + an 1z + an = 0

aplicando la transformada bilineal r se sigue que:

✓ ◆n ✓ ◆n 1
1+r 1+r 1+r
a0 + a1 + · · · + an 1 + an = 0.
1 r 1 r 1 r
n
Simplificando al multiplicar ambos lados de la expresión anterior por (1 r) , se obtiene:

Q(r) = b0 rn + b1 rn 1
+ · · · + bn 1r + bn = 0

y se aplica el criterio de Routh de igual forma como se hace para sistemas continuos.
Ejemplo: Considere la ecuación característica de un sistema de control en tiempo discreto dada como [5]:

z 3 + 5.94z 2 + 7.7z 0.386 = 0.

Sol: Al aplicar la transformada r se obtiene

3.128r3 11.74r2 + 2.344r + 14.17 = 0.

Tabulación de Routh-Hurwitz Considere el siguiente polinomio

a0 s n + a1 s n 1
+ · · · + an 2s
2
+ an 1s + an = 0

35
 sn a0 a2 a4 a6 ···
n 1
s a1 a3 a5 a7 ···
sn 2
b1 b2 b3 b4 ···
sn 3
c1 c2 c3 c4 ···
sn 4
d1 d2 d3 d4 ···
.. .. ..
. . .
s2 e1 e2
s1 f1
s0 g1
donde
a1 a2 a0 a3 b1 a 3 a 1 b2
b1 = c1 =
a1 b1
a1 a4 a0 a5 b1 a 5 a 1 b3 c 1 b2 b1 c 2
b2 = c2 = d1 =
a1 b1 c1
a1 a6 a0 a7 b1 a 7 a 1 b4 c 1 b3 b1 c 3
b3 = c3 = d2 =
a1 b1 c1
.. .. ..
. . .
El criterio de Routh-Hurwitz establece que: Suponiendo a0 > 0: El polinomio P (s) no tiene raíces del
lado derecho del plano complejo s si y solo si todos los pivotes son positivos, es decir b1 > 0, c1 > 0, ..., g1 > 0.
1
Ejemplo: Considere un sistema discreto que tiene una FT G(z) = 2
y un compensador pro-
z (z + z + 1)
porcional con ganancia K en cascada. Determinando el polinomio característico en lazo cerrado simple se tiene
[5]:
z 3 + z 2 + z + K = 0.
Determine el valor para K tal que el sistema discreto sea estable.
Sol: al aplicar la transformada r se obtiene
(1 K) r3 + (1 + 3K) r2 + 3 (1 K) r + 3 + K = 0.

Tabulación de Routh-Hurwitz
r3 1 K 3 (1 K)
r2 1 + 3K 3+K
8K (1 K)
r1 0
1 + 3K
r0 3+K
Las condiciones para que se garantice la estabilidad son:
1 K > 0, 1 + 3K > 0, K > 0, 3 + K > 0.
Por lo tanto
0 < K < 1.

5.2.2. Criterio de Joury

Existen pruebas de estabilidad que se pueden aplicar directamente a las FT en z. Un método tabular sencillo
es el que se conoce como criterio de estabilidad de Joury (1964). Desafortunadamente estas técnicas se vuelen
muy tediosas para sistemas de orden mayor a dos, especialmente cuando se tienen parámetros desconocidos o
por determinar, como el es caso cuando se tiene una ganancia por definir su valor. En la actualidad, no existen
razones suficientes para utilizar estas técnicas si se conocen los coeficientes de la ecuación característica, ya que
siempre se puede utilizar una computadora para determinar las raíces de la ecuación característica e identificar
si hay alguna fuera del circulo unitario.
Previo al análisis del criterio de Joury, se puede partir de las siguientes condiciones necesarias que deben
satisfacerse para que la ecuación característica no tenga raíces sobre o fuera del circulo unitario [5]. Considere
F (z) = an z n + an 1z
n 1
+ · · · + a2 z 2 + a1 z + a0 = 0

36
 F (1) > 0
F ( 1) > 0 si nes entero par
(17)
F ( 1) < 0 si nes entero impar
|a0 | < an .

En otras palabras, si se cumplen las condiciones anteriores, hay que construir el arreglo de Jury para deter-
minar la estabilidad del sistema. Si no se cumplen, se puede dar por hecho que al menos hay una raíz fuera del
circulo unitario [5].
Para un sistema de segundo orden F (z) = a2 z 2 + a1 z + a0 = 0, las condiciones anteriores (17) son necesarias
y suficientes para determinar la estabilidad del sistema, es decir, si se cumplen, el sistema sera estable, de lo
contrario inestable. Ejemplo F (z) = z 2 + z + 0.25 = 0 [5].
Ejemplo: Sea la EC

F (z) = z 3 + z 2 + 0.5z + 0.25 = 0.

Sol: F (1) = 2.75 > 0 y F ( 1) = 0.25 < 0 |a0 | = 0.25 < a3 = 1.

Por lo que se satisfacen las condiciones, sin embargo nada se puede decir sobre la estabilidad del sistema.
Ejemplo: Sea la EC

z 3 + z 2 + 0.5z + 1.25 = 0.

Sol: Como F ( 1) = 0.75 > 0 y como |a0 | = 1.25 no es menor que a3 = 1, entonces la EC tiene al menos
una raíz fuera del circulo unitario.
Tabulación para verificar el Criterio de Joury:
Considere el siguiente polinomio

a0 z n + a1 z n 1
+ · · · + an 2z
2
+ an 1z + an = 0

donde
renglon
1: a0 a1 a2 ··· an 1 an
2: an an 1 an 2 ··· a1 a0
3: b0 b1 ··· bn 2 bn 1
4: bn 1 bn 2 ··· b1 b0
5: c0 ··· cn 3 cn 2
6: cn 2 ··· c1 c0
.. ..
. .
2n + 1 : z0

a0 a0 an an b0 b0 bn 1 b n 1
b0 = c0 =
a0 b0
a1 a0 an 1 an bn 2 b0 b1 bn 1
b1 = cn 2 =
a0 b0
an 1 a0 a1 an
bn 1 =
a0
.. .. ..
. . .

El criterio de Joury establece lo siguiente: Suponiendo que a0 > 0: El polinomio P (z) no tiene raíces
fuera del disco unitario si y solo si todos los pivotes son positivos, es decir b0 > 0, c0 > 0, ..., z0 > 0.
Si ningún pivote es cero, el número de pivotes negativos es igual al número de raíces fuera del disco unitario.
Obsérvese que por cada par de renglones del arreglo anterior solamente se calcula el primero y el siguiente
simplemente se invierte en orden respecto al calculado.

37
 Obsérvese también que si los elementos sombreados (que aquí llamaremos pivotes) en el arreglo de Joury
son cero, los cálculos no se pueden continuar.
Ejemplo: Considere la ecuación característica siguiente

z3 0.1z 2 0.12z 0.4 = 0

y determine por la prueba de Joury si es estable o no.

a0 a1 a2 a3
1 -0.1 -0.12 -0.4
-0.4 -0.12 -0.1 1
0.84 -0.148 -0.16
Sol: Es Estable. Las raíces son: z0 = 0.8281, z1 = 0.364 + 0.5919 j y
-0.16 -0.148 -0.84
0.809 -0.176
-0.176 0.809
0.77
z2 = 0.364 0.5919 j, por tanto estable.
Tarea: Determinar por la prueba de Joury si la EC. (z 2) (z + 0.5) (z 0.1) = 0.
Sol: Es Inestable. Ademas, no se cumple la primera de las condiciones necesarias (F (1) > 0) para no tener
raíces sobre o fuera del circulo unitario.

5.2.3. Lugar Geométrico de las Raíces

La ecuación característica 1 + KG(z)H(z) = 0 debe satisfacer las siguientes dos condiciones:

Magnitud

|KG(z)H(z)| = 1 (18)

Fase o ángulo

\KG(z)H(z) = ±180(2k + 1), k 2 N+ [ 0 = {0, 1, 2, . . .} para K > 0 (19)

Los valores de z que cumplen tanto las condiciones de ángulo como la de magnitud son las raíces de la
ecuación característica, o los polos en lazo cerrado. Ya que comúnmente K varia de cero a infinito, siempre
habrá un valor de K tal que (18) sea cumplida, por lo que la condición que va a restringir un lugar sobre el
LGR es (19). Todo punto que pertenezca al lugar de las raíces debe cumplir la condición de ángulo.

5.2.4. Reglas para la Construcción del Lugar Geométrico de las Raíces de 1 + KG(z)H(z) = 0.

1. Los puntos en K = 0. Los puntos en K = 0 son los polos en G(z)H(z), incluyendo aquellos en z = 1.
2. Los puntos en K = ±1. Los puntos en K = ±1 son los ceros de G(z)H(z) incluyendo aquellos en
z = 1.
3. Numero de los LGR separados. El numero total del lugar geométrico de las raíces es igual al orden
de la ecuación G(z)H(z).
4. Simetría del LGR. Los LGR simétrico respecto al eje real.
5. Asíntotas de LGR cuando z ! 1. Para valores grandes de z, el LR (K > 0) son asintóticos o son
asíntotas con ángulos dados por
2i + 1
✓i = ⇥ 180
|n m|

38
 para el lugar geométrico de las raíces complementario (CRL) (K < 0)

2i
✓i = ⇥ 180
|n m|

en donde i = 0, 1, 2, ..., |n m| 1; n =número de polos finitos de G(z)H(z), m =número de ceros finitos

de G(z)H(z).
6. Intersección de las asíntotas.
a) La intersección de las asíntotas estará solamente sobre el eje real en el plano z.
b) El punto de intersección de las asíntotas es dado por:
P P
partes reales de polos de G(z)H(z) partes reales de ceros de G(z)H(z)
1 =
|n m|

7. Lugar geométrico de las raíces sobre el eje real. Los LR (K > 0) se encuentran en una sección de
eje real sólo si el número total de polos y ceros reales de G(z)H(z) a la derecha de la sección es impar. Si
el número total de polos reales y ceros a la derecha de una sección dada es par, se encuentra CRL (K < 0)
[5]. Al construir los lugares geométricos sobre el eje real, seleccione un punto en éste. Si la cantidad total
de polos y ceros reales a la derecha de este punto de prueba es impar, este punto se encuentra en el lugar
geométrico de las raíces [8].

8. Ángulos de salida. Los ángulos de salida o de llegada de un lugar de raíces desde un polo o un cero de
G(z)H(z) se puede determinar al suponer un punto z1 que esta muy cercano al polo, o cero, y al aplicar
la ecuación:
m
X n
X
\G(z1 )H(z1 ) = \(z1 + zk ) \(z1 + pj )
k=1 j=1
= (2i + 1) ⇥ 180 RL, K>0
= 2i ⇥ 180 CRL K<0

en donde i = 0, ±1, ±2, . . .

9. Intersección del lugar geométrico de las raíces con el círculo unitario. Los puntos de cruce del
lugar geométrico de las raíces sobre el eje imaginario y el valor correspondiente de K se pueden encontrar
al emplear la extensión del criterio de Routh-Hurwitz (modificado).
10. Puntos de ruptura de las trayectorias. Los puntos de ruptura sobre el LGR se determinan al encontrar
las raíces de dK/dz = 0, ó dG(z)H(z)/dz = 0. Estas son condiciones necesarias solamente. Use regla del
cociente4 .
11. Cálculo de los valores de K sobre el LGR. El valor absoluto de K en cualquier punto s1 sobre el
LGR se determina de la ecuación:
1
|K| = .
|G1 (z1 )H1 (z1 )|

Kz 1 e T
Ejemplo: Determine el lugar geométrico de las raíces para el sistema G(z) = considerando
(z 1) (z e T )
un periodo de muestreo de 0.5 segundos [9].
Para el periodo de muestreo mencionado la función de transferencia tiene polos en z = 1 y en z = 0.6065 y
un cero en z = 0.
Para el trazado del lugar de las raíces se ubican primeramente los polos y ceros de la FT y a continuación
los puntos de ruptura de salida y de entrada. Note que esta FT con dos polos y un cero da un LGR circular con
4 d g(z) g 0 (z)h(z) g(z)h0 (z)
=
dz h(z) h2 (z)

39
centro en el origen. El punto de ruptura de salida y el punto de ruptura de entrada se determinan escribiendo
1
la EC. de la formaK = , es decir
G(z)H(z)

(z 1) (z 0.6065)
K=
0.3935z
por tanto
dK z 2 0.6065
=
dz 0.3935z 2
De ahí z 2 = 0.6065, es decir z = 0.7788 y z = 0.7788. Note que la sustitución de los valores anteriores de z en
la función para K dan valores de K = 0.1244 y K = 8.041, respectivamente. El primer valor de z corresponde
al punto de ruptura de salida real mientras que el segundo es el punto de ruptura de entrada real.La siguiente
figura ilustra el LGR

De la ecuación de la ganancia se tienen que al tomar el valor absoluto con el fin de determinar el valor de la
ganancia critica se tiene
(z 1) (z 0.6065)
K=
0.3935z
En vista de que la ganancia critica Kcr corresponde a z = 1, sustituimos ese valor en la ecuación previa
obteniendo
( 2) ( 1.6065)
Kcr = = 8.165.
0.3935( 1)

Tarea: Determine el LGR para la función de transferencia [6]

0.15K (z + 0.7453)
G(z) =
(z 1) (z 0.4119)

donde K es un parámetro variable.

5.2.5. Criterio de Nyquist

Es un método gráfico utilizado para determinar la estabilidad de un sistema en lazo cerrado mediante la
gráfica polar de la función de transferencia en lazo abierto GH(z)[6]. La FT de lazo cerrado de un sistema de
control digital con una entrada y una salida esta dada por [6]

Y (z) G(z)
Gcl (z) = = .
R(z) 1 + G(z)

40
Entonces se puede determinar la estabilidad de un sistema a partir de la ecuación característica
1 + GH(z) = 0.
Semejante a como se realiza para sistemas continuos, el criterio de Nyquist para sistemas de control digital
requiere lo siguiente [6]:

1. Definir la trayectoria de Nyquist en el plano z que contiene el exterior del circulo unitario.
2. Mapear la trayectoria de Nyquist en el plano complejo z sobre el plano G(z) con la función G(z). Con
esto se obtiene la gráfica de Nyquist.
3. La condición de estabilidad del sistema de lazo cerrado se obtiene investigando el comportamiento de la
gráfica de Nyquist de G(z) con respecto al punto crítico ( 1, j0) del plano G(z).

Previo a utilizar con ventaja los métodos desarrollados de la respuesta en frecuencia para sistemas continuos
al análisis y diseño de sistemas de control en tiempo discreto, son necesarias ciertas modificaciones al plano
complejo z [9]. Dado que en el plano z la frecuencia aparece en la forma z = ej!T , donde la exponencial
2 3
T s (T s) (T s)
1+ + + + ···
se puede aproximar como eT s = 2 8 48 y donde al tratar de analizar la respuesta en
2 3
T s (T s) (T s)
1 + + ···
2 8 48
frecuencia del plano z, la simplicidad de las trazas logarítmicas se perderá totalmente. La dificultad, sin embargo,
puede resolverse transformando la función de transferencia de pulso discreta en el plano z en la correspondiente
en el plano w. La transformada llamada comúnmente transformada w, es una transformada bilineal definida
por
1 + (T /2) w
z=
1 (T /2) w
donde T es el periodo de muestreo. Por otro lado, resolviendo para w se tiene que
2 z 1
w= .
T z+1
La transformación puede verse en la siguiente figura.

Ejemplo: Considere la función de transferencia que se muestra en la siguiente figura y un periodo de

muestreo de T = 0.1s. Determine la función de transferencia G(w) [9].

41
 La función de transferencia queda como i considerando que esT representa un retraso de un periodo de
muestreo:
⇢
1 e sT 10
G(z) = Z
s s + 10
⇢
10
= 1 z 1 Z
s (s + 10)
0.6321
= .
z 0.3679
Aplicando la transformación bilineal se tiene que
1 + 0.05w
z=
1 0.05w
por lo que
0.6321
G(z) =
1 + 0.05w
0.3679
1 0.05w
0.6321 (1 0.05 w)
=
0.6321 + 0.06840 w
1 0.05w
= 9.241 .
w + 9.241
Observe que la localización del polo es s = 10 para el sistema continuo, mientras que en el nuevo dominio de
la variable compleja w está en w = 9.241. Similar ocurre en lo que respecta a la ganancia. Ademas aparece un
cero que no estaba presente en el sistema original y mismo que se encuentra ubicado en el semiplano derecho
del plano complejo w (ocasionando inestabilidad a alta ganancia en un sistema de control en lazo cerrado).
Tarea:

Equipo 1: Diagrama de Nyquist, Diagrama de Bode y Margenes de Estabilidad (Incluir ejemplos).

Equipo 2: Especificaciones de Desempeño en Función de los Margenes de Estabilidad.

Equipo 3: Diseño de un Compensador de Atraso Adelanto Usando el Diagrama de Bode (incluir simulación
antes y después del sistema compensado)
Equipo 4: Diseño basado en el LGR.

Diseño Basado en el Lugar Geométrico de las Raíces En el diseño basado en la respuesta en frecuencia
de los sistemas se busca modificar la respuesta en frecuencia tal que se logren determinados margenes de
estabilidad, características de la respuesta transitoria y de estado estable, etc. Aunque el procedimiento basado
en frecuencia tiene metodologías bien establecidas, aun sigue siendo un procedimiento a prueba y error finalmente
dado que comúnmente no se consiguen los resultados deseados en una primera iteración [10].
La compensación de adelanto produce, en esencia, un mejoramiento razonable en el tiempo de la respuesta
transitoria y un cambio pequeño en la precisión en estado estable. Un inconveniente de este esquema es que
puede acentuar los efectos del ruido de alta frecuencia. Por su parte, la compensación de atraso produce un
mejoramiento notable en la precisión en estado estable a costa de aumentar el tiempo de respuesta transitoria.
Suprime los efectos de las señales de ruido a altas frecuencias. La compensación de atraso-adelanto combina las
características de la compensación de adelanto con las de la compensación de atraso. El uso de un compensador
de atraso o de adelanto aumenta el orden del sistema en 1 (a menos que haya una cancelación entre el cero del
compensador y un polo de la función de transferencia en lazo abierto no compensada). El uso de un compensador
de atraso-adelanto eleva el orden del sistema en 2, lo que complica más el análisis y diseño [8] y [9].
El diseño basado en el LGR es una técnica de control que da mejores resultados que los basados en frecuencia.
En términos generales, el LGR es una gráfica que ilustra el movimiento de las raíces de la ecuación característica

42
con respecto a la variación de la ganancia, haciendo evidente la respuesta del sistema en función de la ubicación
de tales raíces. En el procedimiento de diseño consiste en agregar polos y ceros por medio de un compensador
(controlador) tal que obligue a que las raíces de la ecuación característica se desplacen a ubicación mas apropiadas
en el plano z [10].
El procedimiento para el diseño de un compensador de adelanto basado en el LGR es el siguiente [9]:
Considere el sistema mostrado en la siguiente figura

La compensación mediante adelanto es útil cuando el sistema es inestable para todos los valores de ganancia
o es estable, pero tiene características de respuesta transitoria no deseables [9].

1. De las especificaciones de desempeño, determine la posición deseada para los polos dominantes en lazo
cerrado.
2. Calcule la deficiencia angular tal que se tengan los polos en la posición deseada. Este angulo adicional
debe ser proporcionado por el compensador de adelanto tal que el nuevo LGR pase por ahí.
3. El compensador de adelanto tiene la función de transferencia
1 + Tz
C(z) = Kd ↵ , 0 < ↵ < 1.
1 + ↵T z
1
z+
= Kd T
1
z+
↵T
z z0 1 1
= Kd , zp < z0 , z0 = , zp =
z zp T ↵T

4. Si no se especifican las constantes de error estático, determine la localización del polo y del cero del
compensador de adelanto (sobre el eje real y dentro del círculo unitario) de tal forma que el compensador
contribuya con el angulo necesario . Si no hay otro requisitos en el diseño, haga el valor de ↵ tan grande
como sea posible, lo que dará como resultado un valor grande de la constante de velocidad Kv , lo cual
siempre es deseable. (Si se especifica una constante de error estático particular, por lo general siempre es
mas sencillo utilizar el método de diseño basado en frecuencia).
5. Determine la ganancia en lazo abierto del sistema compensado a partir de la condición de magnitud.
6. Finalmente, si no se satisfacen las condiciones de desempeño requeridas, re-ubique tanto el polo como el
cero del compensador a fin de cumplir con los requisitos de desempeño.

En este diseño, el LGR es desplazado hacia la izquierda [10], resultando en un sistema que responde más rápido.
En resumen, el compensador de retraso de fase desplaza muy poco el LGR pero permite tener un ajuste de
ganancia relativamente alto, o si se usa la misma ganancia del sistema de control, el sistema será mas estable
[10].
El procedimiento de diseño del compensador de adelanto es el siguiente [10]:

1. Elija la ubicación de la raíz deseada zb en el plano z.

43
 2. Elija el cero del compensador z0 para cancelar un polo de G(z).
3. Elija ya sea la ganancia K del sistema compensado o bien la ubicación del polo del compensador zp < z0 ,
tal que el compensador sea de adelanto de fase.
4. Determine K o zp a partir de
KD(zb )G(zb ) = 1.

Ejemplo Diseñe un compensador de adelanto de fase en tiempo discreto con función de transferencia [10]
z z0
D(z) = Kd
z zp

para el sistema en tiempo continuo dado como

1
Gp (s) = .
s (s + 1)

Para un periodo de muestreo de T = 0.1 s, la función de transferencia del sistema discreto está dada por
⇢
z 1 1
G(z) = Z
z s2 (s + 1)
0.004837 K (z + 0.9672)
= .
(z 1) (z 0.9048)

Con la intensión de no introducir ganancia adicional en el lazo de control mediante el controlador, debemos
1 zp
hacer que D(1) = 1, por lo que Kd = .
1 z0
De acuerdo al procedimiento descrito previamente, se debe elegir un punto zb de interés, el cual no definimos
ahora pero que consideramos esta a la izquierda del LGR actual (sistema sin compensar). Ahora, se cancela
uno de los polos de G(z), en particular el que está ubicado en 0.9048, por medio de colocar en ese punto el
cero de compensador (notar que hay una cancelación exacta del polo y el cero por lo que podríamos prescindir
de ellos en el diseño de ahora en adelante). Eligiendo el polo del compensador en zp = 0.7, entonces se tiene
que Kd = 3.15. Finalmente eligiendo el punto de ruptura como el valor deseado zb para las raíces del sistema
compensado, se tiene que zb = 0.844, entonces la ganancia K se calcula a partir de la magnitud de la función
de transferencia y el compensador como
1
K =
|D(z)G(z)|
|términos en el denominador {D(z)G(z)}|
=
|términos en el numerador {D(z)G(z)}|
z=zb =0.844

|z 1| |z 0.7|
=
3.15 ⇥ 0.004837 |z + 0.9672|
z=zb =0.844
= 0.814.

Adicionalmente se puede determinar la constante de tiempo del sistema compensado mediante la siguiente
relación z = e T /⌧ . Ahora bien, si se tiene la ubicación del sistema compensado en z = zb = 0.844, entonces
⌧ = 0.59 s. Para fines de comparación, el polo previo a la compensación estaba ubicado en el punto de ruptura
era de z = 0.952, de donde se puede calcular que la constante de tiempo del sistema era de ⌧ = 2.03 s, por lo
que el compensador mejora el tiempo de respuesta del sistema.

5.2.6. Diagramas de Bode

44
5.3. Controlabilidad, Alcanzabilidad y Observabilidad
Controlabilidad: Se dice que un proceso G es controlable si cualquier estado de G puede ser
afectado o controlado en tiempo finito por alguna señal de control u(k) no restringida [6]. Para el
caso de que no se pueda llevar a alguna de las variables de estado a un valor determinado, se dice
que esa variable de estado es no controlable.
Observabilidad: Se dice que un proceso G es observable si cualquier variable de estado de un
proceso afecta alguna de las salidas del proceso. Si alguno de los estados no es observable a partir de
las mediciones de las salidas, se dice que el estado es no observable [6].
Alcanzabilidad: (Es un concepto similar al de controlabilidad) Un sistema es alcanzable si se
puede llevar al estado desde el origen a otro estado cualquiera en tiempo finito [9].

5.3.1. Controlabilidad de sistemas lineales discretos e invariantes en el tiempo

Suponga el siguiente conjunto de ecuación es dinámicas descritas en espacio de estado como

x(k + 1) = A x(k) + B u(k) (20)

y(k) = C x(k) + D u(k)

donde x(k) 2 Rn es un vector de estado de dimensión n, u(k) 2 Rm es un vector de entradas de dimensión

m, y(k) 2 Rp es un vector de salidas de dimensión p; A, B, C, y D son matrices de coeficientes de dimensión
apropiada.

Definición [10]: El sistema (20) es controlable si se garantiza que existe una secuencia de entra-
das u(0), u(1), u(2), . . . , u(N ) que traslade el sistema de cualquier estado inicial x(0) a cualquier
estado final x(N ) en un numero de iteración N finito (tiempo finito).

A continuación se derivan las condiciones tal que el sistema (20) sea controlable. De la ecuación de estado, la
condición inicial dada y de la entrada de control inicial, se tiene que

x(1) = A x(0) + B u(0)

x(2) = A x(1) + B u(1)
= A2 x(0) + A B u(0) + B u(1)
..
.
x(N ) = AN x(0) + AN 1
B u(0) + · · · + A B u(N 2) + B u(N 1)
2 3
u(N 1)
⇥ ⇤6
6 u(N 2) 7 7
= AN x(0) + B AB ··· AN 1
B 6 .. 7
4 . 5
u(0)

De esta forma, conociendo el estado inicial x(0) y el estado final x(N ), las relaciones anteriores se pueden escribir
como 2 3
u(N 1)
⇥ ⇤6
6 u(N 2) 7
7
B AB · · · AN 1 B 6 .. 7 = x(N ) AN x(0)
4 . 5
u(0)
misma que se convierte en un sistema de n ecuaciones lineales en términos de la entrada de control u. Para que
la solución de la ecuación anterior exista, el rango de
⇥ ⇤
B AB · · · AN 1 B

debe ser n, es decir ⇥ ⇤

rank B AB ··· AN 1
B = n.

45
5.3.2. Observabilidad de sistemas lineales discretos e invariantes en el tiempo

Considere el sistema descrito por (20).

Definición [10]: El sistema (20) es observable si se asegura que el estado inicial x(0) puede ser determinado
a partir de las N mediciones y(0), y(1), . . . , y(N 1), para un valor de N finito.
La derivación de la condición para que el sistema sea observable son semejantes a como se realizo en con-
trolabilidad. Por simplicidad, considere que la entrada para todo k es igual a cero, u(k) = 0, entonces se tiene
la salida como

y(0) = C x(0)
y(1) = C x(1)
= C A x(0)
..
.
y(N 1) = C x(N 1)
N 1
= CA x(0)

o bien 2 3 2 3
y(0) C
6 y(1) 7 6 CA 7
6 7 6 7
6 .. 7=6 .. 7 x(0).
4 . 5 4 . 5
y(N 1) C AN 1

Semejante al caso de la controlabilidad, la expresión anterior es un conjunto de n ecuaciones lineales para x(0),
mismas que tiene solución si el rango de la matriz
2 3
C
6 CA 7
6 7
6 .. 7
4 . 5
N 1
CA

es n, es decir 82 39
>
> C >
>
>
<6 7=>
6 CA 7
rank 6 .. 7 = n.
>
> 4 . 5>>
>
: >
;
C AN 1

6. Diseño de Sistemas de Control en Tiempo Discreto

6.1. Introducción

6.2. Controlador PID continuo

Muchas tareas de control son resueltas por controladores tipo PID, mismo que puede ser representado por
la siguiente ecuación [2]  Z
1 t de (t)
u(t) = K e (t) + e (s) ds + Td (21)
Ti dt
donde e(t) = yref (t) y(t) define el error como la diferencia entre el valor de referencia yref (Set Point) y la
salida del proceso o planta a controlar y (la variable medida). K es la ganancia proporcional del controlador,
Ti es la constante de integración o tiempo de reinicio y Td es la constante derivativa.

46
 El mismo controlador puede ser descrito en el dominio de la frecuencia por medio de la siguiente FT:

1
G(s) = KP 1 + + Td s
Ti s
KI KP
= KP + + KD s, KI = , K D = KP T d
s Ti
✓ ◆
T 1
= 1.2 1+ + 0.5L s
L 2L s
✓ ◆2
1
s+
L
= 0.6 T
s
1
por lo que el PID tiene un polo en el origen (un integrador) y un cero doble en s = .
L
Los PIDs fueron implementados inicialmente usando tecnología analógica (válvulas neumáticas, relays y
motores, transistores y circuitos integrados). Actualmente, casi todos los controladores son implementados de
forma digital.
El ultimo termino del controlador PID que corresponde a la parte derivativa, no es fácil de implementar y
de hecho no debería ser implementado dado que este amplificará de forma considerable el ruido existente en la
medición de la variable censada, por lo tanto la ganancia de la parte derivativa debe ser limitada. Esto se puede
realizar por medio de la siguiente aproximación:
sTd
sTd t
1 + sTd /N
la cual realiza la función de derivar a bajas frecuencias y atenúa para las altas frecuencias y por ende atenuando
el ruido que usualmente tiene la característica de tener altas frecuencias. El valor de N esta típicamente en el
rango de 8 a 20 [1]. [2] Dibujar el diagrama de Bode considerando que usualmente Td < 1.
Un método de sintonización del PID propuesto por Ziegler y Nichols, se puede lograr por medio del análisis
de la respuesta del sistema de primer orden ante una entrada tipo escalón (siempre que el sistema no contenga
integradores), cuya respuesta se podría representar en la siguiente figura

K e Ls
La función de transferencia correspondiente seria G(s) = . El análisis de error en estado estable
T s+1
5
puede realizarse haciendo uso del teorema de valor final y de la transferencia del error del sistema en lazo
R(s)
cerrado E(s) = , bajo la consideración de estar analizando el siguiente sistema
1 + C(s)G(s)
5 Teorema del Valor Final:
lı́m f (t) = lı́m sF (s).
t!1 s!0

47
 Analizando el error en estado estable de un sistema sin retardo de primer orden, que tiene una entrada tipo
escalón de amplitud A y un controlador tipo Proporcional Kp
lı́m f (t) = lı́m sE(s)
t!1 s!0
A 1
= lı́m s
s!0 s K
1 + Kp
T s+1
A T s+1
= lı́m s
s!0 s T s + Kp K + 1
A
= .
1 + Kp K

Analizando el error en estado estable de un sistema sin retardo de primer orden, que tiene una entrada tipo
KI
escalón de amplitud A y un controlador tipo Proporcional-Integral Kp + .
s
lı́m f (t) = lı́m sE(s)
t!1 s!0
A 1
= lı́m s ✓ ◆
s!0 s KI K
1 + Kp +
s T s+1
A 1
= lı́m s ✓ ◆
s!0 s Kp s + KI K
1+
s T s+1
A s (T s + 1)
= lı́m s
s!0 s s (T s + 1) + (Kp s + KI ) K
= 0.

La sintonización se basa en la siguiente tabla

Cabe señalar que los valores determinados por la tabla, son solo valores iniciales del PID, posteriormente
hay que realizar los ajustes necesarios a las ganancias a fin de tener el comportamiento deseado. Aunque
dependerá del cada sistema en particular, se recomienda a partir de los valores iniciales proporcionados por la
tabla, aumentar KP para tener un tiempo de subida menor y posterior a eso disminuir KI para disminuir el
sobre impulso. Con respecto a la ganancia KD , modificarla solo en caso de requerir un mejor comportamiento
transitorio o de ser posible eliminarla para evitar amplificación del ruido.

6.3. Controladores PID digitales

El controlador descrito en la ecuación (21) puede ser discretizado usando algunas de las técnicas de aproxi-
mación tanto para las integrales como para las derivadas. debido a la simplicidad del PID, la aproximación por

48
Euler se puede utilizar, resultando en un mecanismo de discretización muy sencillo [2].
Considerando el error e(t) = yref (t) y(t) y dado el compensador del PID como
 Z
1 t de (t)
u(t) = K e (t) + e (s) ds + Td
Ti dt
Z t
K de (t)
= K (yref (t) y(t)) + e (s) ds + KTd .
| {z } Ti | {z dt }
P (t):=
| {z }
I(t):= D(t):=

Analizando por separado cada termino, la parte proporcional se puede describir como

P (kT ) = K (yref (kT ) y(kT )) .

Note que no se requiere ninguna aproximación dado que es estática solamente.

El termino integral
Z t
K
I(t) = e(s)ds
Ti
que se puede escribir como

˙ K
I(t) = e(t)
Ti
que se puede aproximar por Euler como
KT
I [(k + 1) T ] = I(kT ) + (yref (kT ) y(kT )) .
Ti

La parte derivativa dada como

de (t)
D(t) = KTd
dt
dy (t)
= KTd .
dt

s Td
Tomando en consideración la aproximación s Td t para la parte derivativa
1 + s Td /N

1 dy (t)
Dmod (t) = KTd
Td dt
1+s
✓ N ◆
1 dy (t)
= K Td
Td dt
1+s | {z }
N D(t):=
1
= D(t).
Td
1+s
N

1
Como puede notarse, la diferencia entre Dmod (t) y D(t) está dada por el termino , mismo que
Td
1+s
N
eligiendo un valor grande de N , la diferencia se reducirá. Por simplicidad vamos a considerar que la diferencia
entre Dmod (t) y D(t) es pequeña, por lo que podemos escribir

Td dD(t) dy
+ D(t) u KTd
N dt dt

49
es aproximada tomando diferencias hacia atrás, obteniendo
Td D (kT ) D [(k 1) T ] y(kT ) y [(k 1) T ]
+ D(kT ) u KTd
N T T
✓ ◆
Td Td y(kT ) y [(k 1) T ]
+ 1 D(kT ) D [(k 1) T ] u KTd
NT NT T
✓ ◆
Td + N T Td y(kT ) y [(k 1) T ]
D(kT ) D [(k 1) T ] u KTd .
NT NT T

Por lo tanto, multiplicando la expresión anterior por N T y resolviendo para D(kT ) se obtiene
Td KTd N
D(kT ) = D [(k 1) T ] (y(kT ) y [(k 1) T ]) .
Td + N T Td + N T

Considerando todas las partes del PID, la ley de control esta dada como

u(kT ) = P (kT ) + I(kT ) + D(kT ).

6.3.1. Versión PID de Matlab

Otra versión de PID que se analizará enseguida es la versión que viene en el Toolbox de control de
Matlab/Simulink (versión 2012), cuya ecuación es:

1 N
D(s) = P + I . +D
1 s
1+N
s
Usando la aproximación de Euler (diferencia hacia adelante), en donde
z 1
s⇡
T
se obtiene que el PID discreto de Matlab viene dado por

T N
GP ID (z) = P +I +D
z 1 1
1 + NT
z 1
(P + D N ) z 2 + (I T + N P T 2 P 2 D N ) z + P + D N + I N T 2 IT NPT
= .
z 2 + (N T 2) z + (1 N T )

6.4. PID Digital: Takahashi-Chan-Auslander

Existe otra versión de PID digital, y que como para el caso continuo, existe un mecanismo de sintonización
que considera el PID digital que considera el periodo de muestreo. Considere el PID continuo como
 Z
1 t de (t)
u(t) = K e (t) + e (s) ds + Td
Ti dt
mismo que podemos escribir como
 ✓ ◆
1 e0 + e1 ek 1 + ek ek ek 1
u(k) = K e (k) + + ··· + Ts + Td (22)
Ti 2 2 Ts
donde Ts representa el periodo de muestreo. En la ecuación anterior, la integración del error fue remplazada
por una operación de suma del área de trapecios y la derivada mediante una aproximación por diferencias.

50
 KP KI KD
1
P
R (L + Ts )
0.9 1 0.27 Ts
PI ✓ ◆ KI ✓ ◆2
Ts 2 Ts
R L+ R L+
2 2
1.2 1 0.6Ts 0.6
P ID KI ✓ ◆2
R (L + Ts ) 2 Ts R Ts
R L+
2

Cuadro 1: Método de la respuesta transitoria

El algoritmo de Control Digital Directo normalmente aceptado en la practica tienen la diferencia

Uk = u(k) u(k 1) (23)

conocido como algoritmo de velocidad. A partir de tomar un paso anterior para (22), se puede obtener u(k 1).
A partir de (23) se obtiene
 ✓ ◆
1 ek+1 + ek ek 2e(k 1) + ek 2
Uk = K e (k) e(k 1) + Ts + Td .
Ti 2 Ts

Considerando que e(k) = yref (k) y(k) y ademas que la referencia es una constante, por tanto yref (k) =
yref (k 1) = yref (k 2), entonces se obtiene
 ✓ ◆
Ts yk 1 + yk Td
Uk = K y (k 1) y(k) + yref (k) + [2y(k 1) y(k 2) y(k)] .
Ti 2 Ts

1 K Ts K Td
Seleccionando KP = K KI , KI = y KD = se obtiene la expresión
2 Ti Ts

Uk = KP [y (k 1) y(k)] + KI (yref (k) y(k)) + KD [2y(k 1) y(k 2) y(k)] .

La ley de control resultante resulta en

u(k) = Uk + u(k 1). (24)

Nótese que en la última ecuación sólo el término integral incluye la señal de referencia, por
lo tanto en el esquema de velocidad del PID no puede excluirse esta ganancia.

6.4.1. Sintonización de PIDs digitales (Reglas de Takahashi-Chan-Auslander)

Takahashi, Chan y Auslander propusieron en 1970, un conjunto de reglas, que utilizan los dos métodos
propuestos por Ziegler y Nichols para el controlador PID en continuo, a fin de determinar valores aceptables
para KP , KI , KD . La Tabla 1 muestra los valores propuestos para el ajuste de los parámetros del controlador
PID discreto.
Las reglas de ajuste se derivaron para sistemas lineales pero la mayoría de los sistemas reales contienen
algún grado de no linealidad así que es de interés práctico observar qué tan bien trabajan éstas ganancias con
sistemas reales.
K
donde los parámetros6 se determinan de la siguiente curva y considerando que R = .
T
6 Checar algunas restricciones en: Sintonizacion PID Digital Takahasi.pdf, pp. 28

51
 1
KP Td
Ti
Kcr
P
2
1 Kcr Ts
PI 0.45 Kcr KI 0.54
2 Tcr
1 Kcr Ts 3 Kcr Tcr
P ID 0.6 Kcr KI 1.2
2 Tcr 40 Ts
Cuadro 2: Método de la ganancia critica

Si la sintonización se realiza por el método de la ganancia critica se tendrán los siguientes parámetros
donde Kcr es la ganancia requerida para lograr la la oscilación siguiente

y donde Tcr = Pcr es el periodo de la señal.

Trabajo: Utilizando una velocidad de referencia de 100 rad/s para el sistema en lazo cerrado,
un periodo de muestreo de Ts = 0.01 o Ts = 0.001, realice lo siguiente:

1. Implemente en tiempo discreto el compensador de adelanto de fase para el Modelo del Motor de CD dado
como
390s + 6000
D(s) =
s + 30

52
 2. Determine la ecuación en diferencias del PID de Matlab
U (z) Ts N
= D(z) = P + I +D
E(z) z 1 1
1 + N Ts
z 1
donde Ts representa el periodo de muestreo. Use los valores de ganancias proporcional, integral y derivativa
que se determinaron mediante la sintonización por el método del método de Ziegler-Nichols de la respuesta
del modelo continuo. Determine la ecuación en diferencias como u(k) = f (uk , uk 1 , ..., ek , ek 1 , ..., Ts , KP , Ti , Td ).
Use la ecuación en diferencias obtenida para controlar el Motor de CD tanto para el modelo continuo como
para el discreto. Se propone un Ts = 0.01 o Ts = 0.001.
3. Usando las reglas de Takahashi-Chan-Auslander, sintonice el PID siguiente (ecuación del PID discreto)
 Z
1 t de (t)
u(t) = K e (t) + e (s) ds + Td .
Ti dt

Revisar el siguiente material.

http://www.yoquese.com.ar/resources/external/material_evaluacion_y_gestion_de_proyectos/c3.pdf
http://proton.ucting.udg.mx/dpto/tesis/quetzal/CAPITUL4.html

6.4.2. Enredado del integrador (Efecto wind-up) y su eliminación

6.5. Métodos Basados en la Respuesta a la Frecuencia

6.5.1. Compensador de Adelanto de Fase

6.6. Métodos en Espacio de Estados

6.6.1. Retroalimentación de Estado

Los sistemas de control previamente analizados (basados en frecuencia y el LGR), básicamente utilizaban
sólo una señal para fines de control. Parece razonable suponer que si se tiene mas información sobre el sistema
de control (por ejemplo el estado completo) y ésta es utilizada en el diseño, debería entonces tenerse un mejor
desempeño del sistema. Por tanto será necesario tener todo el estado disponible para retroalimentación, es decir,
una retroalimentación de estado (u = K x) [10]. Para la mayoría de los sistemas de control, tener medido el
estado de forma completa es difícil, por lo que se tiene que recurrir al diseño de estimadores a partir de las señales
disponibles. Afortunadamente se pueden realizar los diseños por separado basados en el principio de separación,
de forma tal que se diseñe por un lado el controlador y por otro el estimador (considerando que la convergencia
del observador sea al menos de entre 2 y 5 veces mas rápido que la dinámica del sistema-controlador en lazo
cerrado).

6.6.2. Colocación de Polos por Retroalimentación de Estados

Este procedimiento consiste en ubicar las raíces de la ecuación característica del sistema en un lugar deter-
minado. Para introducir la técnica, considérese el modelo en espacio de estados de un servomotor descrito en
tiempo discreto como [10]

x(k + 1) = A x(k) + B u(k)

 
1 0.0952 0.00484
= x(k) + u(k) (25)
0 0.905 0.0952
⇥ ⇤
y(k) = 1 0 x(k)

donde x1 corresponde a la posición (ángulo) y x2 es la velocidad. Ambas variables de estado pueden ser fácilmente
medibles, por lo que se puede realizar una retro de estado sin mayor problema. Así, se elige una acción de control

53
u(k) como una combinación lineal de los estados dada por
⇥ ⇤
u(k) = K x, K= k1 k2
= k1 x1 k2 x2 .

El sistema en lazo cerrado (por retro de estado) se puede expresar como

 
1 0.0952 0.00484
x(k + 1) = x(k) [k1 x1 (k) + k2 x2 (k)]
0 0.905 0.0952

1 0.00484 k1 0.0952 0.00484 k2
= x(k)
0.0952 k1 0.905 0.0952 k2
| {z }
Ac

donde Ac corresponde la matriz del sistema en lazo cerrado. La ecuación característica a partir de la anterior
representación en espacio de estado se obtiene como det (zI Ac ) = 0, misma que para el ejemplo en cuestión
produce
z 2 + (0.00484k1 + 0.0952k2 1.905) z + 0.00484k1 0.0952k2 + 0.905 = 0.
Ahora bien, se se desea tener ubicados los polos en un determinado lugar del plano z, se esperaría tener un
polinomio característico deseado en términos de los polos deseados, lo cual se puede expresar como

↵d (s) = (z 1 ) (z 2)
= z2 ( 1+ 2) z + 1 2
= 0

por lo que se podrían igualar las los coeficientes de ambos polinomios en z y entonces resolver para k1 y k2
en función de los valores deseados de 1 y 2 . En el caso de sistemas lineales, se puede generalizar lo anterior,
en donde para un sistema de orden n se tendrá un sistema de n ecuaciones lineales con n incógnitas (se debe
resolver para las ganancias ki , i = 1, . . . , n).
Para el caso general, dado el sistema

x(k + 1) = A x(k) + B u(k)

y dado el teorema de Cayley-Hamilton, la matriz polinomial ↵d (A) puede establecerse en términos del los
coeficientes del polinomio característico deseado como

↵c (A) = An + ↵n 1A
n 1
+ · · · + ↵1 A + ↵0 I.

Finalmente, la ubicación de polos puede determinarse por medio de la formula de Ackerman (no se analizará
su deducción aquí, pero ésta es obtenida de considerar que el sistema se puede escribir en la forma canónica
controlable y de ahí partir a su obtención, ver J. E. Ackerman, “Der Entwurf linearer regelungs systems
in Zustansraum”, Regelunstech. Prozess-Datenverarb, Vol. 7, pp. 297-300,1972 . O bien, ver [4]), la
cual es usada para encontrar los valores de K como
⇥ ⇤⇥ ⇤ 1
K= 0 0 ··· 0 1 B AB ··· An 2
B An 1
B ↵c (A).

Los comandos en Matlab para determinar las ganancias por medio de la formula de Ackerman son: acker(A, B, P )
y place(A, B, P ), donde P es un vector con la ubicación deseada de los polos.

Ejemplo Determine la ganancia K tal que la ecuación característica del sistema (25) en lazo cerrado sea

↵c (z) = z 2 1.776 z + 0.819.

54
6.6.3. Observadores

En la sección anterior, se requiere que el estado del sistema sea completamente medible tal que se pueda
realizar la retro de estado u(k) = K x(k). En general, tal medición es impractica y en algunos casos hasta
imposible. En su lugar, esta sección propone una técnica para la estimación de los estados a partir de la
información disponible o medible del sistema y para que posteriormente los estados estimados sean utilizados
para realizar la retro de estados estimados u(k) = K x̂(k). El sistema que estima los estado de otro sistema se
le llama generalmente observador o estimador de estado [10].
Considere el sistema 20 donde A, B y C son matrices de dimensión apropiada y conocidas. Para 20 se
propone un observador como una copia del sistema original en términos del estado estimado y un termino de
corrección en función del error de estimación, con un modelo matemático dado como

x̂(k + 1) = Ax̂(k) + Bu(k) + L (y(k) ŷ(k))

ŷ(k) = C x̂(k)

donde x̂(k) es el estado estimado de x(k).

Para analizar la dinámica del error de observación considere

e(k) = x(k) x̂(k)

Tomando un paso hacia adelante en la dinámica del error se tiene que

e(k + 1) = A x(k) + B u(k) Ax̂(k) Bu(k) L(y ŷ)

= A x(k) A x̂(k) L C x(k) + L C x̂(k)
= A (x(k) x̂(k)) LC (x(k) x̂(k))
= (A LC) (x(k) x̂(k))
= (A LC) e(k).

Por tanto, para que el error de estimación converga a cero, los eigenvalores de A LC deben de estar dentro del
circulo unitario, mismo que se logra seleccionando adecuadamente la ganancia L del observador. La ganancia
del observador se puede determinar utilizando la formula de Ackermann
2 3 1 2 3
C 0
6 CA 7 6 0 7
6 7 6 7
6
L = ↵e (A) 6 .. 7 6 .. 7
. 7 6 . 7.
6 7 6 7
4 C An 2 5 4 0 5
C An 1
1

Ejemplo Determine la ganancia L para el diseño de un observador para el sistema (25), tal que se tenga un
polinomio característico del error de estimación dado como

↵e (z) = z2 1.638 z + 0.671

2
= (z 0.819)

donde tal ubicación de los polos del error de observación producen una respuesta críticamente amortiguada. Use
la fórmula de Ackerman.

Tarea Realizar la simulación del esquema completo controlador-observador para los ejemplos anteriores.

Análisis en lazo Cerrado Sistema-Controlador-Observador Considerando la retro de estado estimada

u(k) = K x̂(k), se analizará el sistema completo, el cual incluye al controlador y observador. Considerando

55
que originalmente la dimensión del estado era n cuando no se tenia al observador, ahora la dimensión será de
2 n dado que el observador es una copia del sistema de control. El sistema en lazo cerrado es

x(k + 1) = A x(k) + B u(k)

x̂(k + 1) = A x̂(k) + B u(k) + L (y(k) ŷ(k))

y dado que u(k) = K x̂(k)

x(k + 1) = A x(k) B K x̂(k)

x̂(k + 1) = A x̂(k) B K x̂(k) + LC (x(k) x̂(k))
= (A BK LC) x̂(k) + L C x(k)

mismo que puede ser reescrito como

  
x(k + 1) A BK x(k)
=
x̂(k + 1) LC A BK LC x̂(k)
| {z }
Af

donde Af es la matriz de estado del sistema en lazo cerrado.

6.6.4. Retroalimentación de la Salida

7. Control Óptimo
7.1. Propiedades

7.2. Control LQR

Considere el siguiente el problema de determinar una acción de control uk tal que el sistema

xk+1 = A xk + B u k

minimice un funcional de costo de la forma

1
X ⇥ ⇤
J= xTk Q xk + uTk R uk .
k=0

Para resolver este problema, bajo el enfoque de Hamilton-Jacobi-Bellman, se plantea el Hamiltoniano como

H(x, u) =

7.3. Filtro de Kalman

7.4. Control LQG

7.5. Identificación por Mínimos Cuadrados

El controlador Gaussiano lineal cuadrático (LQG) es simplemente la combinación del filtro de Kalman
(estimador lineal cuadrático) con el regulador lineal cuadrático (LQR), y esta basado en el principio de separacion

56
8. Proyecto Final
Diseñe e implemente la versión del controlador PID digital de Takahashi-Chan-Auslander para el siguiente
sistema que corresponde a un motor de CD con función de trasferencia:

!(s) 2
= 2
V (s) s + 12 s + 20.02
donde !(s) es la velocidad angular del motor en rad/s y V (s) es el voltaje aplicado al motor en V olts.
Para la realización de la sintonización del PID digital utilice las reglas de Takahashi-Chan-Auslander.

La simulación completa (planta-controlador) debe realizarse en un archivo .m de

Matlab, mediante ecuaciones en diferencias. (NO se aceptará una simulación de Si-
mulink o equivalente)

Reportar:

1. Una introducción y breve descripción de lo que se hará en el proyecto.

2. El procedimiento detallado del diseño y sintonización del PID digital. Considere un periodo de muestreo
de T = 0.01 s.

3. Realice la discretización de la planta continua con el periodo de muestreo previamente referido y conside-
rando que es precedida por un retenedor de orden cero.
4. Simulación del sistema en lazo cerrado sin el PID. Considere una referencia de !ref = 100 rad/s para la
velocidad angular.

5. Simulación del sistema en lazo cerrado con el PID. Considere el mismo valor de referencia del punto
anterior para la velocidad angular.
6. Conclusiones.

57
9. Información Importante
Por definir...

Referencias
[1] Karl J Astrom and Tore Hagglund. PID Cntrollers: Theory, Design and Tuning, 2nd. Ed. USA, 1995.
[2] Karl J Astrom and Bjorn Wittenmark. Computer Controlled Systems: Theory and Design. Prentice-Hall,
Mainland, China, 1997.

[3] Gene F Franklin, J David Powell, and Abbas Emami-Naeini. Control de Sistemas Dinámicos con Retroali-
mentación. Addison Wesley Iberoamericana, Wilmington, Delaware, USA, 1991.
[4] Gene F Franklin, J David Powell, and Michael L Workman. Digital Control of Dynamic Systems. Addison
Wesley Longman, Menlo Park, CA, USA, 1997.
[5] Benjamin C. Kuo. Sistemas de Control Automático. Prentice-Hall Hispanoamericana, Upper Saddle River,
NJ, USA, 1997.
[6] Benjamin C. Kuo. Sistemas de Control Digital. Compania editorial cotinental, Mexico, D.F., 1997.
[7] Frank L Lewis and Vassilis L Syrmos. Optimal Control. John Wiley & Sons, New York, NY, USA, 1995.
[8] K Ogata. Ingenieria de Control Moderno. Prentice-Hall, Madrid, Espana, 2006.

[9] Katsuhiko Ogata. Sistemas de Control en Tiempo Discreto. Prentice-Hall Hispanoamericana, Mexico, 1996.
[10] Charles L Phillips and H Troy Nagle. Digital Control System Analysis and Design. Prentice-Hall Interna-
tional, Inc, Englewood Cliffs, NJ, USA, 1995.

58