AutomaticaII Problemas Parcial1

PROBLEMAS DE

AUTOMATICA II (Primer parcial)

Jose Manuel Díaz Joaquín Aranda

Departamento de Informática y Automática E.T.S.I. Informática U.N.E.D

Problemas de Automática II (Primer parcial). Última reimpresión: septiembre 2008. ISBN: 84-690-0485-9 Copyright  2008 Jose Manuel Díaz y Joaquín Aranda Todos los derechos reservados. Estos apuntes son GRATUITOS y puede ser impresos libremente para fines no comerciales. Sin embargo, la utilización del contenido de estos apuntes en otras publicaciones requiere la autorización de los autores. Departamento de Informática y Automática E.T.S.I Informática. Universidad de Educación a Distancia (UNED). C/ Juan del Rosal nº 16. Madrid 28040 (España)

INDICE ENUNCIADOS............................................................................................................................1 Problemas Tema 1....................................................................................................................1 Problemas Tema 2....................................................................................................................3 Problemas Tema 3....................................................................................................................6 Problemas Tema 4....................................................................................................................8 Problemas Tema 5....................................................................................................................9 SOLUCIONES...........................................................................................................................11 Problema 1.1............................................................................................................................11 Problema 1.2............................................................................................................................12 Problema 1.3............................................................................................................................13 Problema 1.4............................................................................................................................14 Problema 1.5............................................................................................................................16 Problema 2.1............................................................................................................................18 Problema 2.2............................................................................................................................19 Problema 2.3............................................................................................................................22 Problema 2.4............................................................................................................................24 Problema 2.5............................................................................................................................26 Problema 2.6............................................................................................................................27 Problema 3.1............................................................................................................................30 Problema 3.2............................................................................................................................35 Problema 3.3............................................................................................................................41 Problema 3.4............................................................................................................................44 Problema 3.5............................................................................................................................47 Problema 4.1............................................................................................................................49 Problema 4.2............................................................................................................................51 Problema 4.3............................................................................................................................53 Problema 4.4............................................................................................................................54 Problema 5.1............................................................................................................................57 Problema 5.2............................................................................................................................61 Problema 5.3............................................................................................................................64 Problema 5.4............................................................................................................................66 BIBLIOGRAFIA........................................................................................................................74

PROBLEMAS DE AUTOMATICA II (1er parcial)

ENUNCIADOS

ENUNCIADOS PROBLEMAS TEMA 1: Modelos de sistemas continuos 1.1 Una representación en el espacio de estados de un sistema en la forma canónica controlable se obtiene mediante las siguientes ecuaciones:

1   x1  0  x&1   0  x&  = − 0.4 − 1.3· x  + 1·u   2    2  x  y = [0.8 1]· 1   x2 

(1)

El mismo sistema se representa mediante la siguiente ecuación en el espacio de estados, que está en forma canónica observable:

 x&1  0 − 0.4  x1  0.8  x&  = 1 − 1.3 · x  +  1 ·u   2    2  x  y = [0 1]· 1   x2 

(2)

a) Demostrar que la representación en el espacio de estados obtenida mediante las ecuaciones (1) produce un sistema que es de estado controlable pero no observable. b) Demostrar que la representación en el espacio de estados definida mediante las ecuaciones (2) produce un sistema que no es de estado completamente controlable pero si observable. c) Explicar qué provoca las diferencias evidentes en la controlabilidad y la observabilidad del mismo sistema. 1.2 Considere el sistema dado por

 x&1  2 0 0  x1  0 1  x&  = 0 2 0· x  + 1 0· u1   u    2   2   x& 3  0 3 1  x3  0 1  2   x1  1 0 0    y= · x 2  0 1 0  x   3 ¿Es el sistema de estado completamente controlable y completamente observable? ¿Es el sistema de salida completamente controlable? 1.3 Considérese el sistema definido mediante

1   x1  0  x&1   0  x&  = − 2 − 3· x  + 2·u   2    2  x  y = [1 0]· 1   x2  Se desea ubicar los valores característicos del sistema en -3 y -5 usando un control mediante la realimentación del estado u=-K·x. Determinar la matriz de ganancias de realimentación K necesaria y la señal de control u.

1


ENUNCIADOS

1.4 Considérese el sistema definido mediante

1 0   x1  0  x&1   0  x&  =  0 0 1 · x 2  + 0·u  2   x& 3  − 6 − 11 − 6  x3  1  x1  y = [1 0 0]· x 2   x3  Se quiere diseñar un observador de orden completo. Determine la matriz de ganancias del observador Ke usando: a) El método de sustitución directa. b) La fórmula de Ackermann. Supóngase que los valores característicos deseados de la matriz de ganancias del observador son

µ1 = −2 + j·2· 3

µ 2 = −2 − j·2· 3

µ 3 = −5

1.5 Para el sistema definido en el problema 1.4 se supone que la variable de estado x1 es igual a y. En consecuencia x1 es medible y no necesita observarse. Determinar la matriz de ganancias del observador Ke para el observador de orden mínimo. Los valores característicos deseados son

µ1 = −2 + j·2· 3

µ 2 = −2 − j·2· 3

2


ENUNCIADOS

ENUNCIADOS PROBLEMAS TEMA 2: Modelos de sistemas discretos 2.1 Considérese el sistema definido por

G( z ) =

z +1 z + 1.37 z + 0.4 2

Obtener la representación en el espacio de estados para este sistema en las tres configuraciones siguientes: a) Forma canónica controlable b) Forma canónica observable c) Forma canónica diagonal (Jordan) 2.2 Considérese el sistema con doble integrador

T 2 / 2 1 T  x k +1 =  · x + ·u k  k  T 0 1    donde T es el periodo de muestreo. Se pide: a) Determinar la matriz de ganancia de realimentación del estado K tal que la respuesta a una condición arbitraria sea con oscilaciones muertas, es decir, los polos en lazo cerrado deseados deberán estar en el origen de modo que µ1=µ2=0. b) Para el estado inicial x(0)=[1 1]T determinar u(0) y u(T), para T=0.1 seg, T=1 seg y T=10 seg. ¿Que ocurre con las magnitudes de u(0) y u(T) al aumentar T?. 2.3 Considérese la planta definida por

G( z) =

1 z + z + 0.16 2

Utilizando el enfoque de las ecuaciones polinomiales, diseñar un sistema de control para esta planta basado en el diagrama de bloques siguiente R(z )

K0

B( z ) A( z )

U (z )

-

α ( z) F ( z)

Y (z )

β ( z) F ( z)

Suponer que la ecuación característica deseada es H ( z ) = z − 1.2· z + 0.52 y que el polinomio del 2

observador de orden mínimo es F ( z ) = z .

3


ENUNCIADOS

2.4 Sea la planta

1 0  0 0    0 1 ·x k + 0·u x k +1 =  0 − 0.5 − 0.2 1.1 1  x1  y = [0 1 0]· x 2   x3  Usando el enfoque de ecuaciones polinomiales, diseñar un sistema de control para esta planta cuya diagrama de bloques sea R(z )

U (z )

K0

-

B( z ) A( z )

Y (z )

β ( z) α ( z) considérese que

H ( z) = z 3 F ( z) = z 2

2.5 Para la planta del problema 2.4 supuestos los mismos H(z) y F(z) emplear el enfoque polinomial para diseñar un sistema de control para la planta cuyo diagrama de bloques sea R(z )

K0

B( z ) A( z )

U (z )

α ( z) F ( z)

4

-

β ( z) F ( z)

Y (z )


ENUNCIADOS

2.6 Considérese una planta definida por

G( z) =

z + 0.5 z − z + 0.01z + 0.12 3

2

El periodo de muestreo T es de un segundo. Supóngase que en el presente problema es importante tener un error de seguimiento cero a una entrada rampa unitaria. Se quiere diseñar el siguiente sistema de control para la planta R(z )

B( z ) A( z )

U (z )

V (z )

Gm ·H1 ( z )

α ( z) F ( z)

-

β ( z) F ( z)

Se considera que

Gm ( z ) =

0.64·z − 0.512 0.64·z − 0.512 = 3 ( z − 1.2 z + 0.52)·( z − 0.6) z − 1.8·z 2 + 1.24·z − 0.312 2

Se pide: a) Obtener los polinomios α(z) y β(z). b) Obtener las expresiones de Y(z)/V(z) e Y(z)/R(z).

5

Y (z )


ENUNCIADOS

ENUNCIADOS PROBLEMAS TEMA 3: Modelos de perturbaciones 3.1 Calcular la función de covarianza estacionaria del proceso

0  0.4 x k +1 =  · x k + v k − 0.6 0.2 donde v es un proceso de ruido blanco con media cero y varianza

1 0  R1 =   0 2  3.2 Considérese el proceso

− a 0  x k +1 =  · x k + v k  0 − b y k = [1 1]· x k donde v es un proceso de ruido blanco con media cero y matriz de covarianza

σ 12 0  R1 =  2  0 σ2  Demostrar que yk puede ser expresado en la forma

y (k ) = λ ·

q+c e( k ) (q + a)·(q + b)

donde e(k) es ruido blanco con media cero y varianza unidad. Encontrar la relación a partir de la cual es posible determinar λ y c. 3.3 Un proceso estocástico y(k) esta descrito por

x k +1 = a·x k + v k y k = x k + ek donde v y e son procesos de ruido blanco de distribución normal con las propiedades

E[v] = E[e] = 0 Var[v] = 1 Var[e] = r2 r si k = j E[v(k )e( j )] =  12  0 si k ≠ j Demostrar que y(k) puede ser expresada como la salida del siguiente filtro lineal:

6


y (k ) = λ ·

ENUNCIADOS

q−c ε (k ) q−a

| c |< 1

donde ε(k) es ruido blanco con media cero y varianza unidad. Determinar λ y c. 3.4 Determinar la función de covarianza ry(τ) y el espectro φy(ω) del proceso y(k)

y (k ) − 0.7· y (k − 1) = e(k ) − 0.5·e(k − 1) donde e(k) es ruido blanco de varianza unidad. 3.5 Determinar la varianza del proceso estocástico y(k) definido por

y (k ) − 1.5· y (k − 1) + 0.7· y (k − 2) = e(k ) + 0.2·e(k − 1) donde e(k) es ruido blanco de varianza unidad.

7


ENUNCIADOS

ENUNCIADOS PROBLEMAS TEMA 4: Estimación óptima 4.1 Se genera un proceso estocástico mediante

x k +1 = 0.5· x k + v k y k = x k + ek donde v y e son procesos de ruido blanco no correlacionados con covarianzas r1 y r2, respectivamente. Además x(0) tiene una distribución normal con media cero y desviación estándar σ. ¿Cual es la ganancia en estado estacionario?. Calcúlese el polo del filtro estacionario y compárese con el polo del sistema.

4.2 El integrador doble con ruido de proceso puede describirse mediante las ecuaciones:

0  0 .5  1 1 x k +1 =  · x k +  ·u k +    1 0 1 v k  y k = [0 1]x k donde v(k) es una secuencia de variables aleatorias normales de media nula e independientes. Supóngase que x(0) es normal de media nula y matriz de covarianza unidad. Determinar las ecuaciones de la matriz de covarianza del error de estimación y el vector de ganancias del filtro de Kalman.

4.3 Considere el problema anterior pero con salida

y k = [0 1]x k + ek Determinar la covarianza del error de estimación y la ganancia del filtro de Kalman si la varianza de e(k), señal de media nula, es σ2.

4.4 Dado el sistema

0  0 .5  1 1 x k +1 =  · x k +  ·u k +    1 0 1 v k  y k = [0 1]x k donde v(k) es ruido blanco de media nula y desviación estándar 0.1. Supóngase que se conoce x(0) perfectamente. Determinar la estimación de x(k+3), dado y(k) que minimiza el error de predicción.

8


ENUNCIADOS

ENUNCIADOS PROBLEMAS TEMA 5: Identificación de sistemas 5.1 Considérese el sistema:

y (t ) − 0.8· y (t − 1) = u (t − 1) + e(t ) − 0.8·e(t − 1) Se supone que las entradas u(.) y e(.) son ruido blanco de varianza 1. Para este sistema se considera un modelo de la forma:

y (t ) + a· y (t − 1) = b·u (t − 1) a) Obtener por mínimos cuadrados las estimas de los parámetros de este modelo. b) Cuando t→∞ ¿Qué parámetro coincide con el del sistema real y cuál no? ¿Por qué?. c) ¿Se puede conseguir que las estimas de los dos parámetros (a y b) no estén polarizadas? ¿Como?

5.2 Considérese el sistema:

y (k ) + a· y (k − 1) = u (k ) + e(k ) donde u es la variable de control (supóngase que es ruido blanco de varianza λ2), y es la salida y {e(k)} es una secuencia de variables aleatorias gaussianas independientes de varianza σ2. Considérese un controlador autosintonizante basado en el modelo:

y (k ) + β · y (k − 1) = u (k − 1) + e(k ) a) Obtener la estima de mínimos cuadrados de β basada en los datos disponibles hasta el instante k (es decir, y(k), y(k-1),..., y(1), u(k-1),..., u(1)) b) Estudiar si β →a cuando k→∞. Nota: El término e(k) en este modelo está indicando que no existe ruido correlacionado en el modelo. Luego a la hora de realizar los cálculos no hay que tener en cuenta este término.

5.3 Para un sistema descrito por

y (t ) + a· y (t − 1) = b·u (t − 1) + e(t ) con una ley de control dada por

u (t ) = − K · y (t ) ¿Es posible identificar los parámetros de este modelo en lazo cerrado? Justifique razonadamente la respuesta.

9


ENUNCIADOS

5.4 Considérese el sistema:

y (k ) + a· y (k − 1) + b· y (k − 2) = u (k ) + e(k ) donde u es la variable de control, y es la salida y {e(k)} es una secuencia de variables aleatorias gaussianas independientes de varianza σ2. Considérese un modelo de la forma:

y (k ) + θ1 · y (k − 1) + θ 2 · y (k − 2) = u (k − 1) + e(k ) a) Discutir qué tipo de señales de entrada se debe aplicar para realizar la identificación. b) Obtener por mínimos cuadrados las estimas de los parámetros de este modelo. c) Estudiar si θ1→a y θ1→b cuando k→∞. Nota: El término e(k) en este modelo está indicando que no existe ruido correlacionado en el modelo. Luego a la hora de realizar los cálculos no hay que tener en cuenta este término.

10


SOLUCIONES

SOLUCIONES PROBLEMAS TEMA 1: Modelos de sistemas continuos

♦ Solución problema 1.1 a) Para el sistema en la forma (1) la matriz de controlabilidad Mc es:

1  1   0   0 0  0  ·  =  M c = [B | A·B ] =  |  1  − 0.4 − 1.3   1   1 − 1.3  El rango de Mc es 2 ya que |Mc|=-1, luego el sistema es de estado completamente controlable. Por otra parte, la matriz de observabilidad es:

0.8  0 − 0.4   0.8   0.8 − 0.4   M o = C T | AT ·C T =  |  ·  =   1  1 − 1.3   1   1 − 0.5 

[

]

El rango de Mo es 1 ya que |Mo|=0, luego el sistema es no observable. b) Para el sistema en la forma (2) la matriz de controlabilidad Mc es:

0.8  0 − 0.4   0.8   0.8 − 0.4   ·  =  M c = [B | A·B ] =  |   1  1 − 1.3   1   1 − 0.5  El rango de Mc es 1 ya que |Mc|=0, luego el sistema no es de estado completamente controlable. Por otra parte, la matriz de observabilidad es:

1  1   0   0 0  0  M o = C T | AT ·C T =  |  ·  =  1  − 0.4 − 1.3   1   1 − 1.3 

[

]

El rango de Mo es 2 ya que |Mo|=-1, luego el sistema es observable. c) La diferencia aparente en la controlabilidad y la observabilidad del mismo sistema la provoca el hecho de que el sistema original tiene una cancelación de polos y ceros en la función de transferencia asociada al sistema en la forma (1): −1

− 1  0  s G ( s ) = C ·( s·I − A) ·B = [0.8 1]·  ·  = 0.4 s + 1.3 1   s + 1.3 1 0 s + 0.8 1 = 2 ·  = ·[0.8 1]·  s + 1.3·s + 0.4  − 0.4 s  1  ( s + 0.8)·( s + 0.5) −1

11


SOLUCIONES

Nótese que se obtendría la misma función de transferencia si se hubiese considerado el sistema en la forma (2). Si ocurre una cancelación de polos y ceros en la función de transferencia, la controlabilidad y la observabilidad varían, dependiendo de cómo se eligen las variables de estado. Recuérdese que para ser de estado completamente controlable y observable, la función de transferencia no debe tener cancelaciones de polos ni ceros. ♦

♦ Solución problema 1.2 La matriz de controlabilidad del sistema es

0 1 0 2 0 4  M c = B | A·B | A ·B = 1 0 2 0 4 0    0 1 3 1 9 1 

[

2

]

El rango de Mc es 3 ya que

0 1 0  1 0 2  = −3   0 1 3 Luego el sistema es de estado completamente controlable. Por otra parte, la matriz de observabilidad es:

[

( ) ·C

M o = C | A ·C | A T

T

T

2 T

T

]

1 0 2 0 4 0  = 0 1 0 2 0 4    0 0 0 0 0 0 

El rango de Mo es 2, luego el sistema es no completamente observable. Finalmente, para analizar la controlabilidad de la salida hay que estudiar el rango de la siguiente matriz

[

]

0 1 0 2 0 4 0 0 M cy = C·B | C· A·B | C· A 2 ·B | D =   1 0 2 0 4 0 0 0 El rango de Mcy es 2, luego el sistema es de salida completamente controlable. ♦

12


SOLUCIONES

♦ Solución problema 1.3 En primer lugar hay que comprobar que el sistema es de estado completamente controlable, ya que sólo en dicho caso podrá realizarse una ubicación arbitraria de polos. La matriz de controlabilidad para este sistema es:

1   0   0 2  0  0  ·  =  M c = [B | A·B ] =  |  2  − 2 − 3   2    2 − 6  El rango de Mc es 2, ya que |Mc|=-4, luego el sistema es de estado completamente controlable. Por tanto es posible la ubicación arbitraria de polos. Para sistema de ordenes pequeños (n≤ 3) una forma sencilla de calcular K es igualando la ecuación característica en lazo cerrado s·I − ( A − B·K ) ,

 s s·I − ( A − B·K ) =  2 + 2·K1

−1  = s 2 + (3 + 2·k 2 )·s + ( 2 + 2·k1 ) s + 3 + 2·K 2 

con la ecuación característica deseada:

( s + 3)·( s + 5) = s 2 + 8·s + 15 Por tanto,

3 + 2·k 2 = 8 2 + 2·k1 = 15 Resolviendo las ecuaciones anteriores se obtienen los elementos del vector de ganancias K

K = [ k1

k 2 ] = [6.5 2.5]

Asimismo, la señal de control u es:

x  u = − K · x = [6.5 2.5]· 1   x2  ♦

13


SOLUCIONES

♦ Solución problema 1.4 En primer lugar hay que comprobar la observabilidad del sistema, para ello se calcula la matriz de observabilidad:

[

( ) ·C

M o = C | A ·C | A T

T

T

2 T

T

]

1 0 0  = 0 1 0    0 0 1 

Como |Mo|=1 el rango de Mo es 3, por lo tanto el sistema es de estado completamente observable y es posible la determinación de la matriz de ganancias del observador Ke. a) Determinación de Ke mediante sustitución directa Este método consiste en igualar la ecuación característica del observador con la ecuación característica deseada. Para este sistema se tendría:

s·I − A + K e ·C = ( s + µ1 )·( s + µ 2 )·( s + µ 3 ) = ( s + 2 − j·2· 3 )·( s + 2 + j·2· 3 )·( s + 5) = s 3 + 9·s 2 + 36·s + 80

Desarrollando el lado izquierdo de la ecuación anterior se obtiene:

1 0   k e1  s + k e1  s 0 0  0 0 s 0 −  0    0 1 + k e 2 ·[0 0 1] = k e 2       0 0 s  − 6 − 11 − 6  k e3  6 + k e3

−1 s 11

0 −1 = s+6

= s 3 + ( k e1 + 6)·s 2 + (6·k e1 + k e 2 + 11)·s + (11·k e1 + 6·k e 2 + k e3 + 6 ) Luego:

s 3 + ( k e1 + 6)·s 2 + (6·k e1 + k e 2 + 11)·s + (11·k e1 + 6·k e 2 + k e3 + 6) = s 3 + 9·s 2 + 36·s + 80 Igualando los coeficientes de las mismas potencias se obtiene el siguiente sistema de ecuaciones:

k e1 + 6 = 9 6·k e1 + k e 2 + 11 = 36 11·k e1 + 6·k e 2 + k e3 + 6 = 80 Resolviendo el sistema se obtiene

14


SOLUCIONES

 k e1   3  K e = k e 2  =  7       k e3  − 1 a) Determinación de Ke mediante la fórmula de Ackerman La fórmula de Ackerman para la obtención de Ke es:

(

K e = en ·M o−1 ·φ ( AT )

)

T

Para aplicar esta fórmula hay que calcular en primer lugar la inversa de la matriz de observabilidad:

M

−1 o

1 0 0 = 0 1 0   0 0 1 

−1

1 0 0 = 0 1 0   0 0 1 

La ecuación característica deseada para el observador es

∆c ( s ) = ( s + µ1 )·( s + µ 2 )·( s + µ 3 ) = ( s + 2 − j·2· 3 )·( s + 2 + j·2· 3 )·( s + 5) = s 3 + 9·s 2 + 36·s + 80

Por el Teorema de Caley-Hamilton se obtiene

φ ( A) :

φ ( A) = A3 + 9· A2 + 36· A + 80·I Y evaluándola en AT:

φ ( AT ) = (AT ) + 9·(AT ) + 36·(AT ) + 80·I = 3

 −6 = − 11   − 6 74 = 25   3

2

36 − 150  0 − 54 324   0 0 − 216  80 0 0  60 − 239  + 0 − 99 540  + 36 0 − 396  +  0 80 0         25 − 90  9 − 54 225  0 36 − 216   0 0 80  − 18 − 42  − 41 − 95  7 − 1 

Luego sustituyendo en la fórmula de Ackerman y operando se obtiene:

15


SOLUCIONES T

 3 1 0 0 74 − 18 − 42        K e =  [0 0 1]· 0 1 0 · 25 − 41 − 95  =  7         − 1    −  0 0 1 3 7 1     ♦

♦ Solución problema 1.5 El vector de estado y las matrices A y B se pueden reescribir de la siguiente forma:

 x1  − −   x a    x =  x 2  = − −  x  x   3  b  

| 1 0   0 Aab  | − − − − − − − −   Aaa   = − − − − − − A =   0 | 0 1   Abb  |   Aba  | − 11 − 6  − 6

 0   − −   Ba  B =   = − −  0     Bb   1 

Luego se tiene:

x  xb =  2   x3 

x a = x1

Aaa = 0

Aab = [1 0]

0 Aba =    − 6

1  0 Abb =   − 11 − 6

0  Ba = 0 Bb =   1 Para que sea posible la determinación de la matriz de ganancias del observador Ke y poder diseñar el observador de orden mínimo, la condición de observabilidad que se debe cumplir es que:

 A  1 0  O min =  ab  =    Aab · Abb  0 1 sea de rango n-1=3-1=2. Efectivamente esto es así. La ecuación característica deseada para el observador de orden mínimo es:

∆cmin ( s ) = ( s + 10 )·( s + 10 ) = s 2 + 20·s + 100 = 0 Luego

φ ( Abb ) es:

16

(1.1)


SOLUCIONES 2

φ ( Abb ) = Abb

2

1 1  0  0 1 0  89 14  + 100· + 20· Abb + 100·I =  + 20·   =  − 11 − 6 − 11 − 6 0 1  − 154 5 

Luego la matriz de ganancias del observador es:

(

K e = φ ( Abb )· M

)

min −1 o

−1

·e

T n −1

 89 14  1 0 0 14  = ·  ·  =   − 154 5  0 1  1   5  ♦

17


SOLUCIONES

SOLUCIONES PROBLEMAS TEMA 2: Modelos de sistemas discretos

♦ Solución problema 2.1 El sistema tiene la siguiente función de transferencia

G( z ) =

z +1 z + 1.3z + 0.4

G( z ) =

b0 z + b1 z + a1 z + a2

2

Comparando con la expresión:

2

Se obtiene:

b0 = 1, b1 = 1, a1 = 1.3, a2 = 0.4 a) Forma canónica controlable

1  1   0 0   0 0  · xk +  ·uk =  · xk +  ·u k xk +1 =    1  − 0.4 − 1.3 1  − a2 − a1  y k = [b1 b0 ]· xk = [1 1]· xk b) Forma canónica observable

1   0 g  xk +1 =  · x k +  1 ·uk  − a2 − a1  g2  y k = [1 0]· xk Los elementos g1 y g2 provienen del desarrollo en serie de Laurent en el origen de G(z). ∞

G ( z ) = ∑ g k ·z − k = g 0 + g1 ·z −1 + g 2 ·z −2 + ... −1

(2.2)

k =0

En el caso de sistemas discretos este desarrollo se obtiene dividiendo de forma directa el numerador y el denominador de la función de transferencia G(z-1), que en este caso es:

G ( z −1 ) =

z −1 + z −2 1 + 1.3· z −1 + 0.4·z −2

18


SOLUCIONES

Dividiendo

z −1 +

1 + 1.3·z −1 + 0.4·z −2

z −2

− z −1 − 1.3· z −2 − 0.4· z −3 − 0.3·z −2 − 0.4· z −3 + 0.3z −2 + 0.39· z −3 + 0.12·z −4 − 0.01·z −3 + 0.12·z −4

z −1 − 0.3· z −2 +....

: : Se obtiene g1=1 y g2=-0.3, luego

1   1   0 · xk +  xk +1 =  ·uk  − 0.3 − 0.4 − 1.3 y k = [1 0]· xk c) Forma canónica de Jordan

λ 0  1 · x k +1 =  1 x + k  1·u k   0 λ2  y k = [ c1 c2 ]· xk Los coeficientes ci y los autovalores λi i=1,2 se obtienen desarrollando en fracciones simples la función de transferencia G(z):

G( z ) =

c c z +1 5/3 −2/3 = + = 1 + 2 z + 1.3z + 0.4 z + 0.5 z + 0.8 z − λ1 z − λ2 2

Luego

0  1 − 0.5 · x k +  ·uk x k +1 =   − 0.8 1  0 y k = [5 / 3 − 2 / 3]· xk ♦

♦ Solución problema 2.2 a) Para el sistema considerado se tiene que

T 2 / 2 G=   T 

1 T  F=  0 1 

Si se define el vector de ganancias de realimentación como

19


SOLUCIONES

K = [ k1

k2 ]

la ecuación característica en lazo cerrado es:

T2 z·I − ( F ·−G·K ) = z − 1 + 2 ·k1 T ·k1

T2 2 2 − T + ·k 2 = z 2 −  2 − T ·k − T ·k ·z + 1 + T ·k − T ·k = 0 1 2 1 2 2  2 2   z − 1 + T ·k 2

Asimismo, la ecuación característica deseada es

z2 = 0 Igualando los coeficientes de ambas ecuaciones características se obtienen las siguientes ecuaciones:

T2 ·k1 − T ·k 2 = 0 2 T2 1 + ·k1 − T ·k 2 = 0 2 2−

Resolviendo se obtiene:

k1 =

1 3 ; k2 = 2 T 2·T

Luego:

1 K= 2 T

3  2·T 

Aunque no lo pide el problema, a continuación se va a demostrar que la respuesta a las condiciones iniciales es con oscilaciones muertas. Suponga que el estado inicial es:

 x1 (0)  a   x (0 ) = b   2    La ecuación con la realimentación del estado es:

x(( k + 1)·T ) = (F − G·K )· x( k ·T ) sustituyendo valores se obtiene

20


SOLUCIONES

 1 x k T (( 1 )· ) +   2  1  x (( k + 1)·T ) =  1  −  2  T

T  4 · x1 ( k ·T )   1  −   x2 ( k ·T ) 2

Se observa que

 1  x1 (T )   2  x (T ) =  1  −  2  T

T   1 ( 0 ) x     4 · 1 = 2    1 1 −   x2 (0 )  − 2  T

 1  x1 ( 2·T )   2  x ( 2·T ) =  1  −  2  T

T  T   1 a + ·b  a     4 · 4 = 2 1  b   1 1  −  − ·a − ·b  T 2 2 

T  1 T  a + ·b  0    4 · 2 4 = 1  1 1    −  − ·a − ·b 0  2  T 2 

Por lo tanto, queda claro que la respuesta es con oscilaciones muertas. b) En primer lugar se va a determinar u(kT)

1 u( k ·T ) = − K · x( k ·T ) = −  2 T

3  · x( K ·T ) 2·T 

Por lo tanto:

1 u( 0 ) = −  2 T

1 u(T ) = −  2 T

a 3·b 3  a  ·  = − 2 −  2·T  b  T 2·T

T   1 3   2 a + 4 ·b  a b = 2+ ·   1 1 2·T  − ·a − ·b T 2·T  T 2 

En el enunciado se indica que a=1 y b=1, luego

u(0 ) = −

1 3 − 2 T 2·T

u(T ) =

1 1 + 2 T 2·T

Para T=0.1 s se obtiene

u(0) = −115

u(T ) = u(0.1) = 105

Para T=1 s se obtiene

21


SOLUCIONES

u(0) = −2.5

u(T ) = u(1) = 1.5

Para T=10 s se obtiene

u(0) = −0.16

u(T ) = u(10 ) = 0.06

Se observa que para un valor pequeño del periodo de muestreo se hacen grandes u(0) y u(T). Al aumentar el valor del periodo de muestreo T se reducen de forma significativa las magnitudes de u(0) y de u(T). ♦

♦ Solución problema 2.3 Para la planta dada

A( z ) = z 2 + z + 0.16 B( z ) = 1 Por lo tanto a1=1, a2=0.16, b0=0, b1=0, b2=1. Para determinar α(z) y β(z), se debe resolver la siguiente ecuación Diofántica:

α ( z )· A( z ) + β ( z )·B( z ) = F ( z )·H ( z ) = D( z ) Se tiene que

D( z ) = F ( z )·H ( z ) = z·( z 2 − 1.2·z + 0.52 ) = z 3 − 1.2·z 2 + 0.52·z Por lo tanto d0=1, d1=-1.2, d2=0.52, d3=0. Además

α ( z ) = α 0 ·z + α1 β ( z ) = β 0 · z + β1 Para resolver la ecuación Diofántica hay que resolver el sistema de ecuaciones:

E⋅M = D Donde E es la matriz de Sylvester E, que en este caso toma la forma

22


a 2 a 1 E= 1  0

0

b2

a2 a1

b1 b0

1

0

SOLUCIONES

0  0.16 0 b2   1 0.16 = b1   1 1   b0   0 1

1 0 0 1  0 0  0 0

Asimismo el vector D es:

d 3   0  d   0.52   D =  2 =   d1  − 1.2     d 0   1  Luego el sistema a resolver es:

0 0.16  1 0.16  1  1  1  0

1 0 α 1   0  0 1 α 0   0.52  ⋅  =   0 0   β 1   − 1. 2       0 0  β 0   1 

Resolviéndolo se obtiene: α1=-2.2, α0=1, β1=0.352 y β0=2.56 Luego

α ( z ) = z − 2.2 β ( z ) = 2.56·z + 0.352 La función de transferencia en lazo cerrado es:

Y ( z ) K 0 ·B ( z ) 1 = = K0 · 2 R( z ) H ( z) z − 1.2·z + 0.52 Para determinar la constante K0 se requiere que y(∞) en la respuesta al escalón unitario sea igual a uno.

lim y ( k ) = lim(1 − z −1 )·Y ( z ) = lim k →∞

z →1

z →1

K0 K z −1 z · 2 · = 0 =1 z z − 1.2· z + 0.52 z − 1 0.32

Luego

K0 = 8

23


SOLUCIONES

Entonces la función de transferencia en lazo cerrado es:

Y ( z) 8 = 2 R( z ) z − 1.2·z + 0.52 ♦

♦ Solución problema 2.4 En primer lugar hay que obtener la función de transferencia de la planta a partir de la representación mediante variables de estado: −1

0  0   z −1 −1  G p ( z ) = C·(z·I − F ) ·G = [0 1 0]· 0 z − 1  ·0      0.5 0.2 z − 1.1 1  1 1 z B( z ) · z  = 3 = = [0 1 0]· 3 2 2 z − 1.1· z + 0.2· z + 0.5  2  z − 1.1· z + 0.2· z + 0.5 A( z )  z  Luego:

A( z ) = z 3 − 1.1·z 2 + 0.2·z + 0.5 B( z ) = z En consecuencia

a1 = −1.1, a2 = 0.2, a3 = 0.5 b0 = 0, b1 = 0, b2 = 1, b3 = 0 Por otra parte el polinomio característico para el sistema completo es:

D( z ) = F ( z )·H ( z ) = z 3 · z 2 = z 5 Por lo tanto:

d 0 = 1, d1 = 0, d 2 = 0, d 3 = 0, d 4 = 0, d 5 = 0 Para determinar α(z) y β(z), se debe resolver la siguiente ecuación Diofántica:

α ( z )· A( z ) + β ( z )·B( z ) = F ( z )·H ( z ) = D( z ) Donde

24


SOLUCIONES

α ( z ) = α 0 · z 2 + α1 · z + α 2 β ( z ) = β 0 ·z 2 + β1 ·z + β 2 Para resolver la ecuación Diofántica hay que resolver el sistema de ecuaciones:


 a3 a  2 a E= 1 1 0  0

0 a3 a2 a1 1 0

0 0 a3 a2 a1 1

b3 b2 b1 b0 0 0

0 b3 b2 b1 b0 0

0   0.5 0 0   0 0.2 0.5 0   b3  − 1.1 0.2 0.5 = b2   1 − 1.1 0.2 b1   0 − 1.1 1   b0   0 0 1


 d 5  0  d  0   4    d  0  D =  3 =   d 2  0   d1  0      d 0  1  Luego el sistema a resolver es:

0 0  0.5  0.2 0.5 0  − 1.1 0.2 0.5  − 1.1 0.2  1  0 1 − 1.1  0 1  0

0 1 0 0 0 0

0 0 1 0 0 0

0 α 2  0 0 α1  0      0 α 0  0 ⋅  =   1   β 2  0 0  β 1  0      0  β 0  1 

Resolviéndolo se obtiene: α2=-0, α1=1.1, α0=1, β2=-0.55, β1=-0.72, β0=1.01 Luego

α ( z ) = z 2 + 1.1·z β ( z ) = 1.01·z 2 − 0.72· z − 0.55

25

0 1 0 0 0 0

0 0 1 0 0 0

0 0  0  1 0  0


SOLUCIONES

Por lo tanto el controlador diseñado es:

β ( z ) 1.01·z 2 − 0.72·z − 0.55 = α ( z) z 2 + 1.1· z La función de transferencia en lazo cerrado es:

K 0 ·α ( z )·B( z ) K ·α ( z )·β ( z ) K 0 ·( z 2 + 1.1· z )·z K 0 ·( z + 1.1) Y ( z) = = 0 = = R( z ) α ( z )· A( z ) + β ( z )·B( z ) H ( z )·F ( z ) z5 z3 Para determinar la constante K0 se requiere que y(∞) en la respuesta al escalón unitario sea igual a uno.


z →1

z →1

z − 1 K 0 ·( z + 1.1) z · · = 2.1·K 0 = 1 z z3 z −1

Luego

K 0 = 0.4762 Entonces la función de transferencia en lazo cerrado es:

Y ( z ) 0.4762·( z + 1.1) = R( z ) z3 ♦

♦ Solución problema 2.5 Para este problema la función de transferencia de la planta es:

G p ( z) =

z B( z ) = z − 1.1· z + 0.2· z + 0.5 A( z ) 3

2

La ecuación diofántica para este problema es:

α ( z )·( z 3 − 1.1·z 2 + 0.2·z + 0.5) + β ( z )·z = z 5 Esta ecuación es la misma que fue resuelta en el problema 2.4, luego

α ( z ) = z 2 + 1.1·z β ( z ) = 1.01·z 2 − 0.72· z − 0.55

26


SOLUCIONES

La función de transferencia en lazo cerrado es:

Y ( z ) K 0 ·B ( z ) z K = = K 0 · 3 = 20 R( z ) H ( z) z z Para determinar la ganancia K0 se requiere que y(∞) en la respuesta al escalón unitario sea igual a uno.


z →1

z →1

z −1 K0 z · · = K0 = 1 z z2 z −1

Luego

K0 = 1 Entonces la función de transferencia en lazo cerrado es:

Y ( z) 1 = R( z ) z 2 ♦

♦ Solución problema 2.6 a) Para la planta dada

B ( z ) = z + 0.5 A( z ) = z 3 − z 2 + 0.01z + 0.12 Por lo tanto a1=-1, a2=0.01, a3=0.12 b0=0, b1=0, b2=1, b3=0.5. Por otra parte, se observa que no existen factores comunes entre A(z) y B(z). En este problema se puede seleccionar H1(z) como

H 1 ( z ) = z 2 − 1.2 z + 0.52 con el objetivo de cancelar una parte del denominador del modelo Gm(z). Conviene resaltar que esta elección de H1(z) no es única, de hecho existe una infinidad de elecciones posibles. El único requisito necesario es que H1(z) debe ser un polinomio estable de grado n-1 (en este problema n=3). Por otra parte se define

H ( z ) = B( z )·H 1 ( z ) = ( z + 0.5)·( z 2 − 1.2 z + 0.52) = z 3 − 0.7·z 2 − 0.08·z + 0.26 Ahora se escoge a F(z) como

27


SOLUCIONES

F ( z) = z 2 Al igual que ocurría con H1(z), esta elección de F(z) no es única, de hecho existe una infinidad de elecciones posibles. El único requisito necesario es que F(z) debe ser un polinomio estable de grado n-1 (en este problema n=3). Luego D(z) es:

D( z ) = F ( z )·H ( z ) = z 5 − 0.7·z 4 − 0.08·z 3 + 0.26· z 2 Por lo tanto d0=1, d1=-0.7, d2=-0.08, d3=0.26, d4=0, d5=0. Ahora se debe resolver la siguiente ecuación Diofántica:

α ( z )· A( z ) + β ( z )·B( z ) = D( z ) Donde α(z) y β(z) son polinomios en z de segundo orden:

α ( z ) = α 0 · z 2 + α1 · z + α 0 β ( z ) = β 0 ·z 2 + β1 ·z + β 0 Para resolver la ecuación Diofántica hay que resolver el sistema de ecuaciones:


 a3 a  2 a E= 1 1 0  0

0 a3 a2 a1 1 0

0 0 a3 a2 a1 1

b3 b2 b1 b0 0 0

0 b3 b2 b1 b0 0

0  0.12 0 0 0.5 0 0 0   0.01 0.12 0 1 0.5 0     b3   − 1 0.01 0.12 0 1 0.5 =  b2   1 − 1 0.01 0 0 1 b1   0 −1 1 0 0 0    b0   0 0 1 0 0 0


28


SOLUCIONES

d 5   0  d   0   4    d 3   0.26  D= =  d 2  − 0.08  d1   − 0.7      d 0   1  Luego el sistema a resolver es:

0 0 0.5 0 0  α 2   0  0.12  0.01 0.12 0 1 0.5 0  α1   0         − 1 0.01 0.12 0 1 0.5 α 0   0.26   ⋅  =   1 1 0 . 01 0 0 1 −   β 2  − 0.08   0 1 0 0 0   β 1   − 0.7  −1       0 1 0 0 0  β 0   1   0 Resolviéndolo se obtiene: α2=-0.1, α1=0.3, α0=1, β2=0.024, β1=-0.118, β0=0.31 Luego

α ( z ) = z 2 + 0.3·z − 0.1 β ( z ) = 0.31·z 2 − 0.118· z + 0.024 b) La función de transferencia Y(z)/V(z) es:

Y ( z) F ( z )·B( z ) 1 1 = = = 2 V ( z ) F ( z )·B ( z )·H 1 ( z ) H 1 ( z ) z − 1.2 z + 0.52 Puesto que:

V ( z) (0.64· z − 0.512 )·( z 2 − 1.2 z + 0.52 ) 0.64·z − 0.512 = G m ·H 1 ( z ) = = R( z ) ( z 2 − 1.2 z + 0.52 )·( z − 0.6) z − 0.6 Se tiene que la función de transferencia Y(z)/R(z) es

Y ( z ) Y ( z )·V ( z ) 1 0.64·z − 0.512 = = 2 · = Gm R( z ) V ( z )·R( z ) z − 1.2 z + 0.52 z − 0.6 ♦

29


SOLUCIONES

SOLUCIONES PROBLEMAS TEMA 3: Modelos de perturbaciones

♦ Solución problema 3.1 Comparando las ecuaciones del proceso con la expresión general de un proceso lineal discreto estocástico se obtiene:

0  0.4 F=  − 0.6 0.2 En primer lugar se va a calcular la incertidumbre en el estado P(k) que de acuerdo con el Teorema 3.1 viene dada por la ecuación:

P (k + 1) = F·P(k )·F T + R1 Por simplificar la escritura se va a utilizar la siguiente notación:

 p (k ) P (k ) =  11  p 21 (k )

p12 (k )   e(k ) = p 22 (k )  g (k )

f (k ) h(k ) 

Luego

 e(k + 1)  g (k + 1) 

f (k + 1)  0.4 0   e( k ) = ·  h(k + 1)  − 0.6 0.2  g (k )

f (k ) 0.4 − 0.6 1 0 · + 0.2  0 2 h(k )   0

Operando se obtiene:

 e(k + 1)  g (k + 1) 

f (k + 1)  0.16·e(k ) + 1 − 0.24·e(k ) + 0.08· f (k )  =  h(k + 1)  − 0.24·e(k ) + 0.08·g (k ) 0.36·e(k ) − 0.12·g (k ) − 0.12· f (k ) + 0.04·h(k ) + 2

Por lo tanto para obtener las expresiones de e(k), f(k), g(k) y h(k) hay que resolver cuatro ecuaciones en diferencias. Se va a resolver en primer lugar la ecuación en diferencias:

e(k + 1) = 0.16·e(k ) + 1 Tomando la transformada Z se obtiene:

E ( z )·z − e(0)·z = 0.16·E ( z ) +

30

z z −1


SOLUCIONES

Despejando E(z)

E( z) =

e(0)·z z + ( z − 0.16) ( z − 1)·( z − 0.16)

y descomponiendo en fracciones simples el segundo termino se obtiene

1   − 0.16  0.84  e(0)·z E( z) = + + 0.84  ( z − 0.16)  z − 0.16 z − 1    Tomando la transformada Z inversa se obtiene:

1  1  k e(k ) = e(0)·0.16 k −  ·0.16 + 0.84  0.84 

k = 1, 2, 3, ...

Supuesto que e(k0)=e0, entonces:

1  1  k e(k 0 ) = e0 = e(0)·0.16 k0 −  ·0.16 0 + 0.84  0.84  Despejando el valor de e(0)

 1   1  −k e(0) = e0 ·0.16 − k0 +  − ·0.16 0  0.84   0.84  Con lo que:

1  1  k −k e(k ) = e0 ·0.16 k − k0 −  ·0.16 0 + 0.84  0.84 

k = 1, 2, 3, ...

En el estacionario k0→ − ∞ entonces

e( k ) =

1 = 1.1905 0.84

Si se compara esta expresión con la expresión de E(z) se observa que e(k) en el estacionario coincide con el coeficiente de la fracción 1/(z-1). Se va a resolver en segundo lugar la ecuación en diferencias:

f (k + 1) = −0.24·e(k ) + 0.08· f (k )

31


SOLUCIONES

Tomando la transformada Z

F ( z )·z − f (0)·z = −0.24·E ( z ) + 0.08·F ( z ) y despejando F(z) se obtiene:

F ( z) =

f (0)·z 0.24·E ( z ) − z − 0.08 z − 0.08

Sustituyendo el valor de E(z) calculado con anterioridad se obtiene la expresión

α3 α 1 ·z α2 f (0)·z + + + z − 0.08 ( z − 0.08)·( z − 0.16) ( z − 0.08)·( z − 0.16) ( z − 0.08)·( z − 1) α 1 = −0.24·e(0) F ( z) =

0.16 0.84 0.24 α3 = − 0.84

α 2 = 0.24·

Por otra parte como

− c·a − d c·b + d c· z + d = b−a + b−a ( z − a)·( z − b) ( z − a) ( z − b) entonces descomponiendo en fracciones simples a F(z) se obtiene

α2

F ( z) =

α3

α2

α3

α1 2α 1 f (0)·z − + − 0.08 + 0.08 − 0.92 + 0.92 z − 0.08 ( z − 0.08) ( z − 0.16) ( z − 0.08) ( z − 0.16) ( z − 0.08) ( z − 1)

Aplicando la transformada Z inversa se obtiene la siguiente expresión si k=1, 2, 3,...

f (k ) = f (0)·0.08 k − α 1 0.08 k −1 + 2α 1 ·0.16 k −1 −

α2 0.08

0.08 k −1 +

α2 0.08

·0.16 k −1 −

α3 0.92

0.08 k −1 +

α3 0.92

Supuesto que f(k0)=f0 entonces

f (k ) = f (k 0 )·0.08 k − k0 − α 1 0.08 k − k0 −1 + 2α 1 ·0.16 k − k0 −1 −

α2 0.08

0.08 k − k0 −1 +

En el estacionario k0→ − ∞ entonces

f (k ) =

α3 0.92

−

0.24 = −0.3106 0.84·0.92

32

α2 0.08

·0.16 k − k0 −1 −

α3 0.92

0.08 k − k0 −1 +

α3 0.92


SOLUCIONES

Si se compara esta expresión con la expresión de F(z) se observa que f(k) en el estacionario coincide con el coeficiente de la fracción 1/(z-1). Por analogía la ecuación en diferencias

g (k + 1) = −0.24·e(k ) + 0.08· g (k ) tiene la siguiente solución en el estacionario:

g (k ) = −

0.24 = −0.3106 0.84·0.92

Por último queda resolver la ecuación en diferencias

h(k + 1) = 0.36·e(k ) − 0.12·g (k ) − 0.12· f (k ) + 0.04·h(k ) + 2 Tomando la transformada z

H ( z )·z − h(0)·z = 0.36·E ( z ) − 0.12·G ( z ) − 0.12·F ( z ) + 0.04·H ( z ) +

2·z z −1

y despejando H(z) se obtiene:

H ( z) =

h(0)·z 0.36·E ( z ) 0.12·G ( z ) 0.12·F ( z ) 2·z + − − + z − 0.04 z − 0.04 z − 0.04 z − 0.04 ( z − 1)·( z − 0.04)

Sustituyendo las expresiones de E(z), F(z) y G(z) obtenidas con anterioridad se obtiene:

 1   − 0.16  0.84  h(0)·z 0.36  e(0)·z H ( z) = · + 0.84   + + z − 0.04 z − 0.04  ( z − 0.16)  z − 0.16 z − 1       α3 α3  α2 α2    2α 1 α1 0.24 f (0)·z · − − + − 0.08 + 0.08 − 0.92 + 0.92  z − 0.04  z − 0.08 ( z − 0.08) ( z − 0.16) ( z − 0.08) ( z − 0.16) ( z − 0.08) ( z − 1)    +

2·z ( z − 1)·( z − 0.04)

Aunque se podría proceder como en los casos anteriores es más sencillo usar la propiedad observada con anterioridad: La solución de h(k) en el estacionario será el coeficiente de la fracción 1/(z-1). Luego basta con quedarnos con los siguientes términos de la expresión anterior

33


SOLUCIONES

α 0.36 0.24· 3 2·z 0.84 0.92 − + ( z − 0.04)( z − 1) ( z − 0.04)( z − 1) ( z − 1)·( z − 0.04) Descomponiendo en fracciones simples y quedándose únicamente con los términos 1/(z-1) se obtiene

0.24·α 3 2 0.36 0.84·0.96 − 0.92·0.96 + 0.96 ( z − 1) ( z − 1) ( z − 1) Luego

h( k ) =

0.24·α 3 0.36 2 1  0.36 0.24·α 3  − + = − + 2 = · 0.84·0.96 0.92·0.96 0.96 0.96  0.84 0.92 

 1  0.36 0.24 2 = + + 2 = 2.6074 · 0.96  0.84 0.92·0.84  Luego P(k) en el estacionario es:

 1.1905 − 0.3106 P(k ) =   − 0.3106 2.6074  Por otra parte, de acuerdo con el Teorema 3.1, la función de covarianza estacionaria del proceso rx(τ) viene dada por la ecuación:

rx (τ ) = F τ ·P(k ) τ ≥ 0 Operando se obtiene que

τ

0 0.4τ 0   0.4 = F =   τ  − 0.6 0.2 b(τ ) 0.2  τ

donde

0  τ −1 b(τ ) =  − 0 . 6 · 0.4τ − j ·0.2 j ∑  j =0  Finalmente realizando un producto de matrices se obtiene:

34

τ =0 τ >0


SOLUCIONES

 1.1905·0.4τ rx (τ ) =  τ 1.1905·b(τ ) − 0.3106·0.2

− 0.3106·0.4τ

  τ ≥0 − 0.3106·b(τ ) + 2.6074·0.2  τ

♦

♦ Solución problema 3.2 Comparando las ecuaciones del proceso con la expresión general de un proceso lineal discreto estocástico se obtiene:

− a 0  F=   0 − b

C = [1 1]

En primer lugar se va a calcular la incertidumbre en el estado P(k) que de acuerdo con el Teorema 3.1 viene dada por la ecuación:

P (k + 1) = F·P (k )·F T + R1

P(0) = R0

Por simplificar la escritura se va a utilizar la siguiente notación:

 p (k ) P (k ) =  11  p 21 (k )

p12 (k )   e(k ) = p 22 (k )  g (k )

f (k ) h(k ) 

Luego

 e(k + 1)  g (k + 1) 

f (k + 1) − a 0   e(k ) · = h(k + 1)   0 − b   g (k )

f (k ) − a 0  σ 12 · + h(k )   0 − b   0

0  σ 22 

Operando se obtiene:

 e(k + 1)  g (k + 1) 

f (k + 1) a 2 ·e(k ) + σ 12 = h(k + 1)   a·b· g (k )

a·b· f (k )   b ·h(k ) + σ 22  2

Luego para obtener las expresiones de e(k), f(k), g(k) y h(k) hay que resolver cuatro ecuaciones en diferencias. Se va a resolver en primer lugar la ecuación en diferencias:

e(k + 1) = a 2 ·e(k ) + σ 12 Tomando la transformada Z

35


SOLUCIONES

E ( z )·z − e(0)·z = a 2 ·E ( z ) + σ 12 ·

z z −1

y despejando E(z) se obtiene:

E( z) =

e(0)·z z + σ 12 · 2 (z − 1)·( z − a 2 ) (z − a )

Descomponiendo en fracciones simples el segundo termino

 − a2 1  2 e(0)·z 2 1− a 1− a2 + σ · + E( z) =  1 2 z −1 (z − a2 ) z −a 

    

y tomando la transformada Z inversa se obtiene:

e(k ) = e(0)·a

2· k

 σ 12 −  2 1− a

 2·k  σ 12 ·a +  2 1− a 

  

k = 1, 2, 3, ...

Supuesto que e(k0)=e0, entonces:

 σ2 e(k 0 ) = e0 = e(0)·a 2·k0 −  1 2 1− a

 2·k0  σ 12 ·a +  2 1− a 

  

Despejando el valor de e(0)

 σ2 e(0) = e0 ·a − 2·k0 +  1 2 1− a

  σ 12  −  2  1− a

 − 2· k 0 ·a 

Con lo que:

e(k ) = e0 ·a

2·( k − k 0 )

 σ 12 −  2 1− a

 2·( k − k0 )  σ 12 ·a +  2 1− a 

En el estacionario k0→ − ∞ si |a|<1 entonces

e( k ) =

σ 12 1− a2

Por analogía, la ecuación en diferencias

36

  

k = 1, 2, 3, ...


SOLUCIONES

h(k + 1) = b 2 ·h(k ) + σ 22 tiene como solución en el estacionario a:

h( k ) =

σ 22 1− b2

Por otra parte la ecuación en diferencias:

f (k + 1) = a·b· f (k ) tiene como solución

f (k ) = (a·b )

k − k0

En el estacionario k0→ − ∞ si |a·b|<1 entonces f ( k ) = 0 . Por analogía, la ecuación en diferencias

g (k + 1) = a·b·g (k ) tiene como solución en el estacionario a g ( k ) = 0 Luego en el estacionario

 σ 12  2 P (k ) = 1 − a  0 

 0   σ 22  1 − b 2 

Por otra parte, de acuerdo con el Teorema 3.1 la función de covarianza del estado viene dada por:

rxx (k + τ , k ) = Fτ ·P(k ),

τ ≥0

Hay que calcular Fτ,

1 0 − a 0  2 a 2 = F 0 =· F , ·   0 − b, F =· 0 1    0

0  3 − a 3 , F =· b2   0

Luego en el estacionario

37

(− a )τ 0  τ = F ,..., ·   − b3   0

0   (− b )τ 


SOLUCIONES

 (−1)τ ·a τ σ 12  2 rx (τ ) =  1 − a  0 

  τ τ 2  (−1) ·b σ 2  1 − b 2  0

Asimismo la función de covarianza de la salida viene dada por:

 (−1)τ ·a τ σ 12  2 ry (τ ) = C ·rx (τ )·C T = [1 1]· 1 − a  0 

  1 (−1)τ ·a τ σ 12 (−1)τ ·bτ σ 22 + · = (−1)τ ·bτ σ 22  1 1− a2 1− b2 1 − b 2  0

Para considerar además la posibilidad de que τ pueda ser negativo, la expresión anterior se debe escribir en la forma: τ

τ

τ

τ

τ

τ

(−1) ·a σ 12 (−1) ·b σ 22 + ry (τ ) = 1− a2 1− b2 o equivalentemente τ

τ

ry (τ ) = (−1) ·a c1 + (−1) ·b c 2 donde

c1 =

σ 12

c2 =

1− a2

σ 22 1− b2

A continuación hay que calcular la función de densidad espectral de la salida a través de su definición:

Φ y (ω ) =

1 2π

∞

∑ r ( k )e

k = −∞

−ikω

y

Que se puede expresar equivalente en la forma:

Φ y (ω ) =

1 2π

∞ ∞ ∞ −1  1   iβω −ikω  − ikω + + = + + r ( 0 ) r ( k ) e r ( k ) e r ( 0 ) r ( β ) e ry (k )e −ikω  ∑ ∑ ∑ ∑ y y y  y  2π  y β =1 k = −∞ k =1 k =1    

Se van a calcular por separado cada uno de los términos. Para τ=0 se tiene

ry (0) = c1 + c 2

38


SOLUCIONES

Para τ>0 se tiene ∞

∞

k =1

k =1

[

]

(

∞

S1 = ∑ ry (k )e −ikω = ∑ (−1) k ·a k ·c1 + (−1) k ·b k ·c 2 ·e −ikω = c1 ·∑ (−1) k · a·e −iω k =1

)

k

+ c 2 ·∑ (−1) k ·(b·e −iω ) ∞

k

k =1

Considerando de forma independiente los términos pares de los impares es posible expresar la ecuación anterior en la forma:

(

∞

S1 = c1 ·∑ a·e −iω

)

2k

k =1

∞

+ c 2 ·∑ (b·e −iω ) − c1 ·∑ ·(a·e −iω ) ∞

∞

2k

k =1

∞

− c 2 ·∑ (b·e −iω ) ∞

2 k −1

k =1

∞

2 k −1

k =1

∞

S1 = c1 ·∑ (a 2 ·e −i 2ω ) + c2 ·∑ (b 2 ·e −i 2ω ) − c1 ·a −1 ·e iω ∑ ·(a 2 ·e −i 2ω ) − c2 ·b −1 ·e iω ∑ (b 2 ·e −i 2ω ) k

k =1

k

k =1

k

k =1

∞

k

k =1

∞

S1 = ( c1 − c1 ·a −1 ·e iω )∑ (a 2 ·e −i 2ω ) + (c2 − c2 ·b −1 ·e iω )·∑ (b 2 ·e −i 2ω ) k

k =1

k

k =1

Se tienen series de progresiones geométricas que poseen una expresión analítica para su suma, con lo que

S1 = (c1 − c1 ·a −1 ·e iω )

a 2 ·e −i 2ω b 2 ·e − i 2ω −1 iω ( c c · b · e )· + − 2 2 1 − a 2 ·e −i 2ω 1 − b 2 ·e −i 2ω

a 2 c1 − c1 ·a·e iω b 2 c 2 − c 2 ·b·e iω S1 = + e i 2ω − a 2 e i 2ω − b 2 Para τ<0 se tiene ∞

∞

β =1

β =1

[

]

S 2 = ∑ ry ( β )e iβω = ∑ (−1) β ·a β ·c1 + (−1) β ·b β ·c 2 ·e iβω siguiendo un procedimiento de cálculo similar al descrito con anterioridad se puede demostrar que:

a 2 c1 − c1 ·a·e − iω b 2 c 2 − c 2 ·b·e − iω S2 = + e − i 2ω − a 2 e − i 2ω − b 2 Luego la función de densidad espectral de la salida es:

a 2 c1 − c1 ·a·e −iω b 2 c 2 − c 2 ·b·e − iω a 2 c1 − c1 ·a·e iω b 2 c 2 − c 2 ·b·e iω  1  Φ y (ω ) = + + + c1 + c 2 +  2π  e − i 2ω − a 2 e − i 2ω − b 2 e i 2ω − a 2 e i 2ω − b 2 

39


SOLUCIONES

O equivalentemente

F ( z) =

1 2π

 a 2 c1 − c1 ·a· z −1 b 2 c 2 − c 2 ·b· z −1 a 2 c1 − c1 ·a·z b 2 c 2 − c 2 ·b·z  c c + + + + +  1  2 z −2 − a 2 z −2 − b 2 z2 − a2 z 2 − b2  

Sacando factor común a c1 y c2,

F ( z) =

1 2π

  a (a − z −1 ) a ( a − z )   b(b − z −1 ) b(b − z )      c · 1 c · + + + + 2  1 2 1 + −2 −2 2 2  2 2 2  z a z a z b z b − − − −     

desarrollando los denominadores,

F ( z) =

1 2π

   a (a − z −1 ) a(a − z )  b(b − z −1 ) b(b − z )      + + + + + c · 1 c · 1  1 2  −1 −1 −1 −1   ( z − b)( z + b) ( z − b)( z + b)    ( z − a )( z + a ) ( z − a )( z + a ) 

simplificando,

F ( z) =

F ( z) =

1 2π

 1   a a  b b   + c 2 ·1 − −1  − − c1 ·1 − −1 + 2π   ( z + a ) ( z + a )  ( z b ) + ( z b )  

  ( z −1 + a )·( z + a ) − a·( z + a ) − a·( z −1 + a )   ( z −1 + b)·( z + b) − b·( z + b) − b·( z −1 + b)     + c · c ·  1  2  ( z −1 + a )·( z + a ) ( z −1 + b)·( z + b)      F ( z) =

1 2π

 c1 ·(1 − a 2 ) c 2 ·(1 − b 2 )  · +  −1  −1  ( z + a )·( z + a ) ( z + b)·( z + b) 

y dejando un denominador común se llega a la siguiente expresión para la función de densidad espectral:

1  c1 ·(1 − a 2 )·( z −1 + b)·( z + b) + c 2 ·(1 − b 2 )·( z −1 + a )·( z + a )  F ( z) = ·  2π  ( z −1 + a )·( z + a )·( z −1 + b)·( z + b)  Se sabe del enunciado que la expresión del filtro H(z) que teniendo como entrada ruido blanco genera la salida y cuya densidad espectral es F(z), es:

H ( z) = λ·

z+c ( z + a)·( z + b)

Por el teorema de factorización espectral se sabe que

40


SOLUCIONES

2π ·F ( z ) = H ( z )·H ( z −1 ) y sustituyendo el valor de H(z) se tiene:

2π ·F ( z ) =

λ2 ·( z + c )·(z −1 + c ) ( z + a)·( z + b)·( z −1 + a )·( z −1 + b)

Comparando esta expresión con la función de densidad espectral obtenida con anterioridad se obtiene la ecuación

λ2 ·(z + c )·(z −1 + c ) = c1 ·(1 − a 2 )·( z −1 + b)·( z + b) + c 2 ·(1 − b 2 )·( z −1 + a)·( z + a) que es equivalente a:

λ2 ·( z + c )·(z −1 + c ) = σ 12 ·( z + b)·( z −1 + b) + σ 22 ·( z + a)·( z −1 + a) Esta es la relación pedida a partir de la cual es posible determinar λ y c ♦

♦ Solución problema 3.3 El proceso estocástico de este problema es similar al del Ejemplo 3.2 de los apuntes, por lo tanto la función de covarianza estacionaria del estado es

rx (τ ) =

τ

r1 ·a 1− a2

y el valor medio del estado en el estacionario es 0. En este problema r1=1 Se va a calcular la función de covarianza estacionaria de la salida y:

ry (τ ) = E[ y k +τ · y k ] = E[( x k +τ + ek +τ )·( x k + ek )] = E[ x k +τ ·x k ] + E[ x k +τ ·ek ] + E[ek +τ ·x k ] + E[ek +τ ·ek ] =

= rx (τ ) + E[a·x k +τ −1 ·ek ] + E[v k +τ −1 ·ek ] + E[ek +τ ·x k ] + re (τ ) Como el estado es independiente del ruido blanco que afecta a la entrada entonces la expresión anterior se reduce a

ry (τ ) = rx (τ ) + E[v k +τ −1 ·e k ] + re (τ ) que de acuerdo con los datos del enunciado es equivalente a

41


SOLUCIONES

 1  + r2 2  1− a  a + r12 ry (τ ) =  1− a2  τ  a  1− a2

si

τ =0

si

τ =1

si

τ >1

A continuación hay que calcular la función de densidad espectral de la salida a través de su definición:

Φ y (ω ) =

1 2π

∞

∑ r ( k )e

k = −∞

−ikω

y

Que se puede expresar equivalente en la forma: ∞ 1  −2  −ikω − iω iω r ( k ) e r ( 1 ) e r ( 0 ) r ( 1 ) e ry (k )e −ikω  = + − + + ∑ ∑ y y y y  2π  k = −∞ k =2  ∞ ∞ 1  − iω −ikω  iβω iω =  ∑ ry (− β )e + ry (−1)e + ry (0) + ry (1)e ∑ ry (k )e  = 2π  β = 2 k =2  k β ∞ 1  a  a  iω  1   a  − iω ∞ a −ikω  iβω e r e r r e e · · = + + + + + + + ∑ = ∑ 12 2 12 2 2 2 2 2 1 − a  1 − a  1 − a  2π  β = 2 1 − a k =2 1 − a  ∞ ∞ k 1  iω β = + α 2 ·e iω + α 1 + α 2 ·e −iω + b∑ a·e −iω  b ∑ a·e 2π  β = 2 k =2 

Φ y (ω ) =

(

)

(

)

donde

1 1− a2 1 + r2 α1 = 1− a2 a + r12 α2 = 1− a2 b=

Puesto que se trata de series de progresiones geométricas, entonces existe una expresión analitica para su suma:

(

)

(

)

2 2 1  a·e iω a·e −iω  iω − iω Φ y (ω ) = + α 2 ·e + α 1 + α 2 ·e − b − b  2π  1 − a·e iω 1 − a·e −iω 

La expresión anterior es equivalente a:

42


F ( z) =

SOLUCIONES

b·a 2  1  − b·a 2 −1 α α α z z + + + − · · 2 1 2   ( z − a )·z  2π  ( z −1 − a )·z −1

Si se define

α 3 = −b·a 2 entonces:

F ( z) =

 α3 α3 1  + + α 2 ·( z + z −1 ) + α 1   −1 −1 2π  ( z − a )·z ( z − a )·z 

Dejando un denominador común se obtiene la siguiente expresión:

F ( z) =

1 α 3 ·z·( z − a) + α 3 ·z −1 ·( z −1 − a ) + α 2 ·( z + z −1 )( z −1 − a )·( z − a ) + α 1 ·( z −1 − a )·( z − a)    2π  ( z −1 − a)·( z − a) 

Se sabe del enunciado que la expresión del filtro H(z) que teniendo como entrada ruido blanco genera la salida y cuya densidad espectral es F(z), es:

H ( z) = λ·

z−c z−a

Por el teorema de factorización espectral se sabe que

2π ·F ( z ) = H ( z )·H ( z −1 ) y sustituyendo el valor de H(z) se tiene:

2π ·F ( z ) =

λ2 ·( z − c )·(z −1 − c ) ( z − a)·( z −1 − a)

Comparando esta expresión con la expresión de F(z) calculada con anterioridad se obtiene la siguiente ecuación:

λ2 ·(z − c )·(z −1 − c ) = α 3 ·z·( z − a) + α 3 ·z −1 ·( z −1 − a) + α 2 ·( z + z −1 )( z −1 − a)·( z − a) + α 1 ·( z −1 − a)·( z − a) que es equivalente a

λ2 ·(1 + c 2 ) − λ2 ·c·(z −1 + z ) = [α 3 − a·α 3 ]·( z 2 + z −2 ) + [−a·α 3 + α 2 ·(1 + a) − a·α 1 ]·( z + z −1 ) + [−2·a·α 2 + α 1 ·(1 + a 2 )]

43


SOLUCIONES

Igualando los coeficientes asociados a las mismas potencias se obtiene el siguiente par de ecuaciones:

λ2 ·(1 + c 2 ) = α 1 ·(1 + a 2 ) − 2·a·α 2 λ2 ·c = a·α 3 − α 2 ·(1 + a) + a·α 1 Que se pueden expresar equivalentemente en la forma:

λ2 ·(1 + c 2 ) = p1 λ2 ·c = p 2 donde

p1 = α 1 ·(1 + a 2 ) − 2·a·α 2 p 2 = a·α 3 − α 2 ·(1 + a) + a·α 1 Resolviendo este par de ecuaciones se obtiene como soluciones:

c=

p1 p12 ± −1 2· p 2 4· p 22

λ=±

p2 c

De estas cuatro soluciones solamente serían válidas aquellas que fuesen reales y cumplan que |c|<1. ♦

♦ Solución problema 3.4 En primer lugar hay que determinar la expresión del filtro H(z) que excitado con la entrada e genera la salida y, para ello hay que aplicar la transformada Z sobre la ecuación en diferencias del enunciado:

Y ( z ) − 0.7·z −1 ·Y ( z ) = E ( z ) − 0.5·z −1 ·E ( z ) con lo que

H ( z) =

Y ( z ) z − 0.5 = E ( z ) z − 0.7

En segundo lugar hay que determinar la expresión temporal de este filtro, es decir h(k), para ello hay que aplicar la transformada Z inversa a la expresión H(z), con lo que se obtiene:

1 h( k ) =  k −1 0.2 ⋅ 0.7

44

si k = 0 si k = 1, 2, 3,...


SOLUCIONES

Por otra en el dominio temporal la relación entre la entrada y la salida viene dada por la expresión: ∞

y (k ) = ∑ h(n)·e(k − n) n =0

Tomando valores medios

∞  m y (k ) = E [ y (k )] = E ∑ h(n)·e(k − n) =  n =0  ∞

∞

n =0

n =0

= ∑ h(n)·E [e(k − n)] = ∑ h(n)·me (k − n) Supuesto que el ruido blanco de la entrada tiene valor medio 0, entonces el valor medio de la salida también sería 0. Por lo tanto, de la definición de covarianza se obtiene:

[

]

ry (τ ) = E y (k + τ )· y T (k ) = T  ∞  ∞   = E  ∑ h(n)·[e(k + τ − n)]· ∑ h(l )·[e(k − l )]  =   l =0    n =0 ∞

[

∞

]

= ∑∑ h(n)·E e(k + τ − n)·e T (k − l ) ·h T (l ) = n =0 l =0 ∞

∞

= ∑∑ h(n)·re (τ + l − n)·h T (l ) n =0 l =0

Desarrollando el sumatorio más interno se obtiene: ∞

∞

∞

n =0

n =0

n =0

ry (τ ) = ∑ h(n)·re (τ + 0 − n)·h T (0) + ∑ h(n)·re (τ + 1 − n)·h T (1) + ∑ h(n)·re (τ + 2 − n)·h T (2) + ...

Puesto que la entrada de la entrada es ruido blanco su función de covarianza es no nula únicamente cuando se evalúa en 0, y en dicho caso de acuerdo el enunciado vale 1, por ello dentro de cada sumatorio solamente existe un término no nulo:

ry (τ ) = h(τ )·h T (0) + h(τ + 1)·h T (1) + h(τ + 2)·h T (2) + .... La expresión anterior es equivalente a: ∞

ry (τ ) = ∑ h(τ + k )·h T (k ) k =0

Y como se trata de un sistema escalar:

45


SOLUCIONES ∞

ry (τ ) = ∑ h(τ + k )·h(k ) k =0

Si τ=0, entonces ∞

∞

∞

k =0

k =1

k =1

ry (0) = ∑ h 2 (k ) = 1 + ∑ (0.2 ⋅ 0.7 k −1 ) 2 = 1 + 0.2 2 ⋅ 0.7 − 2 ∑ (0.7 2 ) k Esta sumatoria es una progresión geométrica de razón 0.72, que tiene una expresión analitica para su suma, por lo tanto:

ry (0) = 1 + 0.2 2 ⋅ 0.7 − 2

0 .7 2 0.2 2 55 1 = + = 2 2 51 1 − 0.7 1 − 0.7

Si τ>0, entonces ∞

ry (τ ) = h(τ ) + ∑ h(τ + k )·h (k ) = 0.2 ⋅ 0.7 T

τ −1

k =1

∞

(

= 0.2 ⋅ 0.7τ −1 + 0.2 2 ·0.7 − 2 ·0.7τ ∑ 0.7 2

∞

+ ∑ 0.2 ⋅ 0.7τ + k −1 ·0.2 ⋅ 0.7 k −1 = k =1

)

k

= 0.2 ⋅ 0.7τ −1 + 0.2 2 ·0.7 − 2 ·0.7τ

k =1

= 0.2 ⋅ 0.7

τ −1

0.7 2 = 1 − 0.7 2

τ

 0.2 0.2 ·0.7 0.2 2  130 + = + ·0.7τ = ·0.7τ 2 2  357 1 − 0.7  0.7 1 − 0.7  2

Luego la función de covarianza de la salida es:

 55  ry (τ ) =  51 130 τ  ·0.7  357

τ =0 τ ≠0

La función de densidad espectral se puede calcular mediante la expresión dada por el teorema 3.3

φ y (ω ) = H (e iω )·φu (ω )·H T (e − iω ) Sustituyendo valores se obtiene:

( (

)( )(

1 e iω − 0.5 e − iω − 0.5 φ y (ω ) = · iω · 2π e − 0.7 e −iω − 0.7 o equivalentemente

46

) )


F ( z) =

SOLUCIONES

( (

1 ( z − 0.5) z −1 − 0.5 · · 2π ( z − 0.7 ) z −1 − 0.7

) ) ♦

♦ Solución problema 3.5 En primer lugar hay que determinar la expresión del filtro H(z) que excitado con la entrada e genera la salida y, para ello hay que aplicar la transformada Z sobre la ecuación en diferencias del enunciado:

Y ( z ) − 1.5·z −1 ·Y ( z ) + 0.7·z −2 ·Y ( z ) = E ( z ) + 0.2·z −1 ·E ( z ) con lo que

H ( z) =

Y ( z ) B( z ) z 2 + 0.2· z = = 2 E ( z ) A( z ) z − 1.5· z 1 + 0.7

Si se supone que:

B( z ) = b0 ·z 2 + b1 ·z + b2 A( z ) = a 0 ·z 2 + a1 ·z + a 2 entonces b0=1, b1=0.2, b2=0, a0=1, a1=-1.5 y a2=0.7. De acuerdo con el Teorema 3.5 la densidad espectral de la señal y esta dada por

φ (ω ) =

1 B( z )·B( z −1 ) 2·π A( z )· A( z −1 )

También se sigue de dicho Teorema 3.5 que la varianza de la señal y está dada por la integral compleja

[ ]

Ey

2

π

π

B ( z )·B ( z −1 ) dz 1 1 = ∫ φ (ω )·dω = · ∫ φ (ω )·e −iω d (e iω ) = · i −π 2·π ·i ∫ A( z )· A( z −1 ) z −π

Puesto que el orden del denominador es 2, se puede demostrar que en dicho caso la varianza viene dada por la expresión:

[ ]

E y2 =

B0 ·a 0 ·e1 − B1 ·a 0 ·a1 + B2 ·(a12 − a 2 ·e1 ) a 0 ·[(a 02 − a 22 )·e1 − (a 0 ·a1 − a1 ·a 2 )·a 1 ]

donde

47


SOLUCIONES

B0 = b02 + b12 + b22 B1 = 2·(b0 ·b1 + b1 ·b2 ) B2 = 2·b0 ·b2 e1 = a 0 + a 2 Sustituyendo los valores de los coeficientes de A y B en las expresiones anteriores se obtiene

B0 = 1.04, B1 = 0.4, B2 = 0, e1 = 1.7 y

[ ]

E y2 =

2.368 37 = = 12.3333... 0.192 3 ♦

48


SOLUCIONES

SOLUCIONES PROBLEMAS TEMA 4: Estimación óptima

♦ Solución problema 4.1 Considerando que la ecuación general de un proceso discreto es:

x k +1 = Φ· x k + Γ·u k + v k y k = C · x k + ek se obtiene que para el proceso del problema Φ=0.5, Γ=0 y C=1. La ganancia de Kalman en el estado estacionario es:

(

K e = P·C T · R2 + C ·P·C T

)

−1

donde

P = Φ·[ P − P·C T ·( R2 + C ·P·C T ) −1 ·C ·P ]·Φ T + R1 En este caso se trata de un sistema escalar luego:

Ke =

P·C r2 + C 2 ·P

 Φ 2 ·r2 ·P P 2 ·C 2  P = Φ 2 · P − + r = ·+ r1  1 (r2 + C 2 ·P)  r2 + C 2 ·P  Expresando la ecuación de la covarianza del error de estimación como una ecuación de segundo orden:

P·r2 + C 2 ·P 2 = Φ 2 ·r2 ·P + r1 ·r2 + r1 ·C 2 ·P C 2 ·P 2 + (r2 − Φ 2 ·r2 − r1 ·C 2 ) P − r1 ·r2 = 0 C 2 ·P 2 + [(1 − Φ 2 )r2 − r1 ·C 2 ]P − r1 ·r2 = 0 y sustituyendo valores se obtiene:

P 2 + [0.75·r2 − r1 ]·P − r1 ·r2 = 0

49


SOLUCIONES

La solución positiva de esta ecuación es:

P=

r1 − 0.75·r2 + (0.75·r2 − r1 ) 2 + 4·r1 ·r2 2

Luego la ganancia de Kalman en el estado estacionario es:

r1 − 0.75·r2 + (0.75·r2 − r1 ) 2 + 4·r1 ·r2 P K = = r2 +·P r1 + 1.25·r2 + (0.75·r2 − r1 ) 2 + 4·r1 ·r2 e

Por otra se pide calcular el polo del filtro estacionario y compararlo con el del sistema. De forma general, el sistema real viene dado por:

x k +1 = Φ·x k + Γ·u k + v k Por su parte el estimador del sistema implementado por el filtro de Kalman viene dado por la ecuación:

xˆ (k + 1 | k + 1) = xˆ (k + 1 | k ) + K ke+1 ·[ y k +1 − C · xˆ (k + 1 | k )] donde:

xˆ (k + 1 | k ) = Φ·xˆ (k | k ) + Γ·u k Sustituyendo xˆ ( k + 1 | k ) en la expresión de xˆ ( k + 1 | k + 1)

xˆ (k + 1 | k + 1) = Φ·xˆ (k | k ) + Γ·u k + K ke+1 ·[ y k +1 − C ·Φ· xˆ (k | k ) − C ·Γ·u k ] y reordenando términos se obtiene:

xˆ (k + 1 | k + 1) = [1 − K ke+1 ·C ]·Φ· xˆ (k | k ) + [1 − K ke+1 ·C ]·Γ·u k + K ke+1 · y k +1 Que en el estacionario toma la forma:

xˆ (k + 1 | k + 1) = [1 − K e ·C ]·Φ· xˆ (k | k ) + [1 − K e ·C ]·Γ·u k + K ke+1 · y k +1 Para este problema, el sistema real es:

x k +1 = 0.5·x k + v k

50


SOLUCIONES

cuyo polo es z=0.5. Mientras que el estimador en el estado estacionario es

xˆ (k + 1 | k + 1) = [1 − K e ]·0.5· xˆ (k | k ) + K e · y k +1 cuyo polo es

z = [1 − K e ]·0.5 Es decir el polo del estimador en el estado estacionario modifica el polo del sistema real con un factor

[1 − K e ] . ♦

♦ Solución problema 4.2 Comparando la expresión del enunciado con la ecuación general de un proceso discreto

x k +1 = Φ· x k + Γ·u k + wk y k = C ·x k + ek donde

wk = Ω·v k se obtiene que

1 1 0.5 0  Φ= Γ=  Ω=   0 1 1 1  C = [0 1] Por otra parte la varianza de wk es:

0 0  0  T T T R1 = E[ wk ·wk ] = E[Ω·v k ·v k ·Ω T ] = Ω·E[v k ·v k ]·Ω T = r1 ·Ω·Ω T = r1 · ·[0 1] =   1  0 r1  Por otra parte como ek es 0 para este sistema, entonces R2=0; Las ecuaciones de la matriz de covarianza del error de estimación y del vector de ganancias del filtro de Kalman son:

P (k + 1 | k ) = Φ·P(k | k )·Φ T + R1

51


SOLUCIONES

(

K ke+1 = P(k + 1 | k )·C T · R2 + C ·P(k + 1 | k )·C T

)

−1

P (k + 1 | k + 1) = P (k + 1 | k ) − K ke+1 ·C ·P(k + 1 | k ) Sustituyendo valores se obtiene:

1 1 1 0 0 0  P (k + 1 | k ) =  ·P(k | k )·   + 0 1 1 1 0 r1 

K

e k +1

0   0   = P(k + 1 | k )· · [0 1]·P(k + 1 | k )·   1  1 

−1

P (k + 1 | k + 1) = P(k + 1 | k ) − K ke+1 ·[0 1]·P(k + 1 | k ) Considerando la siguiente notación

p12 (k )  p 22 (k )

 p (k ) P (k | k ) =  11  p 21 (k )

 p (k + 1) P (k + 1 | k + 1) =  11  p 21 (k + 1)

p12 (k + 1)  p 22 (k + 1)

Entonces operando se obtienen las ecuaciones buscadas

p (k ) + p 21 (k ) + p12 (k ) + p 22 (k ) p12 (k ) + p 22 (k ) P (k + 1 | k ) =  11 p 21 (k ) + p 22 (k ) p 22 (k ) + r1  

K

e k +1

 p (k) + p 22 (k)  =  12 1  p 22 (k) + r1 

p11 (k + 1) = p11 (k ) + p 21 (k ) + p12 (k ) + p 22 (k ) −

T

(p12 (k ) + p 22 (k )) ⋅ (p 21 (k ) + p 22 (k )) (p 22 (k ) + r1 )

p12 (k + 1) = p 21 (k + 1) = p 22 (k + 1) = 0 ♦

52


SOLUCIONES

♦ Solución problema 4.3 Este problema es análogo al anterior pero ahora R2=σ2, en consecuencia las ecuaciones de la matriz de covarianza del error de estimación y del vector de ganancias del filtro de Kalman ahora son:

1 1 1 0 0 0  ·P(k | k )· P (k + 1 | k ) =    + 0 1 1 1 0 r1 

K

e k +1

0   0   = P(k + 1 | k )· · σ 2 + [0 1]·P(k + 1 | k )·   1  1  

−1

P (k + 1 | k + 1) = P(k + 1 | k ) − K ke+1 ·[0 1]·P(k + 1 | k ) Luego operando se obtienen las ecuaciones buscadas:

p (k ) + p 21 (k ) + p12 (k ) + p 22 (k ) p12 (k ) + p 22 (k ) P (k + 1 | k ) =  11 p 21 (k ) + p22(k ) p 22 (k ) + r1  

K

e k +1

 p (k) + p 22 (k) =  12 2  p 22 (k) + r1 + σ

p 22 (k) + r1   p 22 (k) + r1 + σ 2 

T

p11 (k + 1) = p11 (k ) + p 21 (k ) + p12 (k ) + p 22 (k ) −

(p12 (k ) + p 22 (k )) ⋅ (p 21 (k ) + p 22 (k )) (p 22 (k ) + r1 + σ 2 )

p12 (k + 1) = p12 (k ) + p 22 (k ) −

(p12 (k ) + p 22 (k )) ⋅ (p 22 (k ) + r1 ) (p 22 (k ) + r1 + σ 2 )

p 21 (k + 1) = p 21 (k ) + p 22 (k ) −

(p 21 (k ) + p 22 (k )) ⋅ (p 22 (k ) + r1 ) (p 22 (k ) + r1 + σ 2 )

p 22 (k + 1) = p 22 (k ) + r1 −

(p 22 (k ) + r1 ) 2 (p 22 (k ) + r1 + σ 2 ) ♦

53


SOLUCIONES

♦ Solución problema 4.4 Puesto que el estado inicial x(0) se conoce perfectamente, eso significa que la matriz de covarianza del error de estimación en el instante inicial k=0 es:

0 0  P (0 | 0) =   0 0  De acuerdo con las ecuaciones del filtro de Kalman:

0 0  P (1 | 0) = Φ·P(0 | 0)·Φ Y + R1 = R1 =   0 r1  0   0   K = P(1 | 0)· · [0 1]·P(1 | 0)·   1   1  

−1

e 1

0  =  1 

0 0  P (1 | 1) = P(1 | 0) − K 1e ·[0 1]·P(1 | 0) =   0 0  Con lo que se deduce que en cualquier instante k se tendrá:

0 0  P (k | k ) =   0 0  0  K ke+1 = K e =   1  Por otra parte el estimador del sistema implementado por el filtro de Kalman viene dado por la ecuación:

xˆ (k + 1 | k + 1) = xˆ (k + 1 | k ) + K ke+1 ·[ y k +1 − C · xˆ (k + 1 | k )] donde:

xˆ (k + 1 | k ) = Φ·xˆ (k | k ) + Γ·u k Sustituyendo xˆ ( k + 1 | k ) en la expresión de xˆ ( k + 1 | k + 1) :

xˆ (k + 1 | k + 1) = Φ·xˆ (k | k ) + Γ·u k + K ke+1 ·[ y k +1 − C ·Φ· xˆ (k | k ) − C ·Γ·u k ]

54


SOLUCIONES

y reordenando términos se obtiene:

xˆ (k + 1 | k + 1) = [ I − K ke+1 ·C ]·Φ· xˆ (k | k ) + [ I − K ke+1 ·C ]·Γ·u k + K ke+1 · y k +1 Que se puede escribir de forma equivalente como:

xˆ (k + 1 | k + 1) = A· xˆ (k | k ) + B·u k + K ke+1 · y k +1 donde

A = [ I − K ke+1 ·C ]·Φ B = [ I − K ke+1 ·C ]·Γ Para este sistema la ganancia de Kalman es constante, luego

xˆ (k + 1 | k + 1) = A· xˆ (k | k ) + B·u k + K e · y k +1 En el instante k+2 el estimador es:

xˆ (k + 2 | k + 2) = A· xˆ (k + 1 | k + 1) + B·u k +1 + K e · y k + 2 Sustituyendo en la ecuación anterior la expresión del estimador en k+1 y operando se obtiene:

xˆ (k + 2 | k + 2) = A 2 · xˆ (k | k ) + A·B·u k + A·K e · y k +1 + B·u k +1 + K e · y k + 2 En el instante k+3 el estimador es:

xˆ (k + 3 | k + 3) = A·xˆ (k + 2 | k + 2) + B·u k + 2 + K e · y k +3 Sustituyendo en la ecuación anterior la expresión del estimador en k+2 y operando se obtiene:

xˆ (k + 3 | k + 3) = A 3 ·xˆ (k | k ) + A 2 ·B·u k + A 2 ·K e · y k +1 + A·B·u k +1 + A·K e · y k + 2 + B·u k + 2 + K e · y k + 3 Reordenando términos se obtiene:

xˆ (k + 3 | k + 3) = A 3 ·xˆ (k | k ) + [ A 2 ·B·u k + A·B·u k +1 + B·u k + 2 ] + [ A 2 ·K e · y k +1 + A·K e · y k + 2 + K e · y k +3 ]

Por tanto se ha conseguido expresar el estimador en el instante k+3 en función de: el estimador en el instante k, las señales de control en los instantes k, k+1 y k+2, y las salidas medidas en los instantes k+1, k+2 y k+3.

55


SOLUCIONES

Para el sistema del problema las matrices A y B toman los siguientes valores:

0.5 1 1  A= B=   0 0 0  Además se puede comprobar que

An = A

n = 1, 2, 3,...

que

A·B = B y que

1  A·K e =   = D 0  Luego el estimador en el instante k+3 para este sistema toma la expresión:

xˆ (k + 3 | k + 3) = A· xˆ (k | k ) + [ B·u k + B·u k +1 + B·u k + 2 ] + [ D· y k +1 + D· y k + 2 + K e · y k +3 ] 0  1 0.5 1 1 xˆ (k + 3 | k + 3) =  ·xˆ (k | k ) +  ·(u k + u k +1 + u k + 2 ) +  ·( y k +1 + y k + 2 ) +  · y k +3  1 0  0 0 0  ♦

56


SOLUCIONES

SOLUCIONES PROBLEMAS TEMA 5: Identificación

♦ Solución problema 5.1 a) El proceso a identificar es:

y (t ) = 0.8· y (t − 1) + u(t − 1) + e(t ) − 0.8·e(t − 1) El modelo que se considera para la identificación es:

y (t ) = − a· y (t − 1) + b·u(t − 1) En este caso:

a 

φ T (t ) = [ − y (t − 1) u(t − 1)], θ =   b  

Para construir la ecuación normal hay que hacer las siguientes multiplicaciones de matrices

− y (1 − 1)  · − y (1 − 1) · − y (i − 1) · − y (t − 1)  T Φ ·Φ =  · − y (i − 1)  u (1 − 1) · u (i − 1) · u (t − 1)   ·  − y (t − 1)

u (1 − 1)  ·  u (i − 1)   ·  u (t − 1) tx1

 y (1)  ·   − y (1 − 1) · − y (i − 1) · − y (t − 1)    Φ T ·Y =  · y ( i )    u (1 − 1) · u (i − 1) · u (t − 1)    ·   y (t ) que se pueden expresar en la forma: t t   2 ( 1 ) y i y (i − 1)·u (i − 1) − − ∑ ∑  i =1 Φ T ·Φ =  t i =1 · t 2 − y (i − 1)·u (i − 1)  u (i − 1) ∑  ∑  i =1 i =1

57

tx1


SOLUCIONES

 t  − ∑ y (i )· y (i − 1) Φ T ·Y =  ti =1 ·  y (i )·u (i − 1)   ∑  i =1 Luego la ecuación normal es: t t    t  2 − − − − y ( i 1 ) y ( i 1 )· u ( i 1 ) ∑ ∑  a  − ∑ y (i )· y (i − 1)  i =1  t i =1 ·  =  ti =1  t b  2 − y (i − 1)·u(i − 1)   u (i − 1) y (i )·u(i − 1)  ∑ ∑  ∑    i =1 i =1   i =1

Resolviendo el sistema de ecuaciones se obtienen las estimas de mínimos cuadrados de los parámetros a y b

 t  t   t  t  −  ∑ u 2 (i − 1) · ∑ y (i )· y (i − 1)  +  ∑ y (i − 1)·u(i − 1) · ∑ y (i )·u(i − 1)    i =1   i =1   i =1  a =  i =1 2 t t t       ∑ y 2 (i − 1) · ∑ u 2 (i − 1)  −  ∑ y (i − 1)·u(i − 1)   i =1   i =1   i =1   t  t   t  t  −  ∑ y (i − 1)·u(i − 1) · ∑ y (i )· y (i − 1)  +  ∑ y 2 (i − 1) · ∑ y (i )·u(i − 1)    i =1   i =1   i =1  b =  i =1 2 t t t       ∑ y 2 (i − 1) · ∑ u 2 (i − 1)  −  ∑ y (i − 1)·u(i − 1)   i =1   i =1   i =1  b) Multiplicando ambos miembros de la ecuación por 1/t y haciendo la aproximación:

1 t 2 1 t ·∑ x (i − 1) ≈ ·∑ x 2 (i ) t i =1 t i =1 la ecuación normal se puede expresar en la forma:

1 t  1 t 2   1 t  · y ( i ) · y ( i )· u ( i ) − − ·∑ y (i)· y (i − 1) ∑  t∑   a   t i =1 t i =1 i =1   ·  =   t t t 1 1 1 b  2 − · y (i)·u (i)   ·∑ u (i) ·∑ y (i )·u (i − 1)   t ∑    t i =1 i =1   t i =1 Si t→ ∞ , entonces se obtiene

 E[ y 2 (t )] − E[ y (t )·u (t )] a  − E[ y (t )· y (t − 1)]   ·  =   E[u 2 (t )]  b   E[ y (t )·u (t − 1)]  − E[ y (t )·u (t )]

58

(1)


SOLUCIONES

Hay que calcular los valores esperados de diferentes magnitudes. Como la entrada u y el ruido e son independientes se cumple:

E[e(t )·e(t − 1)] = 0 E[u (t )·u (t − 1)] = 0 E[e(t )·u (t )] = 0 Por otra parte

E[ e(t ) 2 ] = 1 E [u(t ) 2 ] = 1 Además como el valor actual de e o de u es independiente del valor pasado de la salida y, se cumple:

E[ y (t − 1)·u (t )] = 0 E[ y (t − 1)·e(t )] = 0 Teniendo en cuenta toda esta información, en primer lugar, se va a calcular E[ y (t )·u (t )] :

E[ y (t )·u(t )] = E [( 0.8· y (t − 1) + u(t − 1) + e(t ) − 0.8·e(t − 1))·u(t )] = = E[0.8· y (t − 1)·u(t )] + E [u(t − 1)·u(t )] + E[ e(t )·u(t )] + E[ −0.8·e(t − 1)·u(t )] = = 0+0+0+0 = 0 En segundo lugar se va a calcular E[ y (t )·u (t − 1)]

E[ y (t )·u(t − 1)] = E[( 0.8 y (t − 1) + u(t − 1) + e(t ) − 0.8e(t − 1))·u(t − 1)] = = E[0.8· y (t − 1)·u(t − 1)] + E[u(t − 1)·u(t − 1)] + E[ e(t )·u(t − 1)] + E[ −0.8·e(t − 1)·u(t − 1)] = = 0 +1 + 0 + 0 = 1 2

En tercer lugar se va a calcular E[ y (t ) ]

E[ y (t ) 2 ] = E[( 0.8· y (t − 1) + u(t − 1) + e(t ) − 0.8·e(t − 1)) 2 ] desarrollando el cuadrado, tomando el valor esperado, y considerando únicamente los términos no nulos se obtiene la siguiente expresión

E[ y (t ) 2 ] = 0.64·E[ y 2 (t − 1)] − 1.28·E[ e(t − 1)· y (t − 1)] + E[u 2 (t − 1)] + + E[ e 2 (t )] + 0.64·E[ e 2 (t − 1)] Se observa que para poder obtener E[ y (t ) ] se requiere calcular E[e(t − 1)· y (t − 1)] : 2

59

(2)


SOLUCIONES

E[ y (t − 1)·e(t − 1)] = E [( 0.8· y (t − 2) + u(t − 2) + e(t − 1) − 0.8·e(t − 2))·e(t − 1)] = = E[0.8· y (t − 2)·e(t − 1)] + E[u(t − 2)·e(t − 1)] + E[ e(t − 1)·e(t − 1)] + E [ −0.8·e(t − 2)·e(t − 1)] = = 0 + 0 +1 + 0 = 1 Asimismo si se consideran las siguientes aproximaciones:

E[ y (t − 1) 2 ] ≈ E [ y (t ) 2 ] E[u(t − 1) 2 ] ≈ E [u(t ) 2 ] = 1 E[ e(t − 1) 2 ] ≈ E[ e(t ) 2 ] = 1 entonces la ecuación (2) toma la forma

E[ y(t ) 2 ] = 0.64·E[ y (t ) 2 ] − 1.28 + 1 + 1 + 0.64 = 0.64·E[ y (t ) 2 ] + 1.36 2

Luego, despejando E[ y (t ) ] se obtiene:

E [ y 2 (t )] =

1.36 = 3.78 0.36

Finalmente se va a calcular E[ y (t )· y (t − 1)]

E[ y (t )· y (t − 1)] = E[( 0.8· y (t − 1) + u(t − 1) + e(t ) − 0.8e(t − 1))· y (t − 1)] = = E[0.8 y (t − 1)· y (t − 1)] + E [u(t − 1)· y (t − 1)] + E[ e(t )· y (t − 1)] + E[ −0.8·e(t − 1)· y (t − 1)] = = 0.8·E [ y (t ) 2 ] + 0 + 0 − 0.8 = = 0.8·3.78 − 0.8 = 2.22

Por lo tanto, la ecuación normal toma la forma:

3.78 0 a  2.22  · =  0 1  b   1   Resolviéndola se obtienen las siguientes estimas:

a = 0.59 b =1 Solamente el parámetro b=1 coincide con el parámetro real b0=1. Por su parte el parámetro a=0.59 difiere del parámetro real a0=0.8 debido a la existencia de ruido correlacionado que hace que la estima esté polarizada. c) Si es posible, para ello se deben utilizar otros métodos de estimación como el método de los

60


SOLUCIONES

mínimos cuadrados extendidos, el algoritmo de máxima verosimilitud o el método de la variable instrumental. ♦

♦ Solución problema 5.2 a) El proceso a identificar es:

y (t ) = −ay (t − 1) + u (t ) + e(t ) donde e(t) es ruido gaussiano independiente. El modelo que se considera para la identificación es:

yˆ (t ) = − β y (t − 1) + u (t − 1) En este caso:

β  φ T (t ) = [− y (t − 1) u (t − 1)], θ =   1  

Para construir la ecuación normal hay que hacer las siguientes multiplicaciones de matrices


u (1 − 1)  ·  u (i − 1)   ·  u (t − 1) tx1

 y (1)  ·   − y (1 − 1) · − y (i − 1) · − y (t − 1)    · Φ T ·Y =  y i ( )    u (1 − 1) · u (i − 1) · u (t − 1)    ·   y (t ) tx1 que se pueden expresar en la forma: t t   2 − − ( 1 ) y i y (i − 1)·u (i − 1) ∑ ∑  i =1 Φ T ·Φ =  t i =1 · t 2 − y (i − 1)·u (i − 1)  u (i − 1) ∑  ∑  i =1 i =1

61


SOLUCIONES

 t  − ∑ y (i )· y (i − 1) Φ T ·Y =  ti =1 ·  y (i )·u (i − 1)   ∑  i =1 Luego la ecuación normal es: t t    t  2 − − − − y ( i 1 ) y ( i 1 )· u ( i 1 ) ∑ ∑    β  − ∑ y (i )· y (i − 1) i =1  t i =1 ·  =  ti =1  t 1  2 − y (i − 1)·u (i − 1)   u (i − 1) y (i )·u (i − 1)  ∑ ∑  ∑    i =1 i =1   i =1

Despejando β se obtiene: t

β=

t

∑ y(i − 1)·u (i − 1) − ∑ y(i)· y(i − 1) i =1

i =1

t

∑y

2

(1)

(i − 1)

i =1

b) Multiplicando el númerador y el denominador de (1) por 1/t y haciendo la aproximación:

1 t 2 1 t ·∑ x (i − 1) ≈ ·∑ x 2 (i ) t i =1 t i =1 la ecuación (1) si t→ ∞ se puede expresar en la forma:

β=

E[ y (t − 1)·u(t − 1)] − E[ y (t )· y (t − 1)] E[ y 2 (t )]

(2)


E[e(t )·e(t − 1)] = 0 E[u (t )·u (t − 1)] = 0 E[e(t )·u (t )] = 0 Por otra parte

E[e(t ) 2 ] = σ 2 E[u (t ) 2 ] = λ2 Además como el valor actual de e o de u es independiente del valor pasado de la salida y, se cumple:

62


SOLUCIONES

E[ y (t − 1)·u (t )] = 0 E[ y (t − 1)·e(t )] = 0 Teniendo en cuenta toda esta información, en primer lugar, se va a calcular E[ y (t )·u (t )] :

E[ y (t )·u (t )] = E[(− a· y (t − 1) + u (t ) + e(t ))·u (t )] = = E[−a· y (t − 1)·u (t )] + E[u (t ) 2 ] + E[e(t )·u (t )] = = 0 + E[u (t ) 2 ] + 0 = λ2 2

En segundo lugar se va a calcular E[ y (t ) ]

E[ y (t ) 2 ] = E[(−a· y (t − 1) + u (t ) + e(t )) 2 ] desarrollando el cuadrado, tomando el valor esperado, y considerando únicamente los términos no nulos se obtiene la siguiente expresión

E[ y (t ) 2 ] = a 2 ·E[ y 2 (t − 1)] +·E[u 2 (t )] + E[e 2 (t )] Asimismo si se consideran las siguientes aproximaciones:

E[ y (t − 1) 2 ] ≈ E[ y (t ) 2 ] entonces la ecuación (3) toma la forma

E[ y (t ) 2 ] = a 2 ·E[ y 2 (t )] +·E[u 2 (t )] + E[e 2 (t )] 2

Luego, despejando E[ y (t ) ] se obtiene:

E[ y (t )] = 2

λ2 + σ 2 1− a2

Finalmente se va a calcular E[ y (t )· y (t − 1)]

E[ y (t )· y (t − 1)] = E[(−a· y (t − 1) + u (t ) + e(t ))· y (t − 1)] = = E[−a· y (t − 1)· y (t − 1)] + E[u (t )· y (t − 1)] + E[e(t )· y (t − 1)] = = − a·E[ y (t ) 2 ] − a·λ2 + −a·σ 2 = 1− a2 Por lo tanto, la ecuación (2) toma la forma:

63

(3)


SOLUCIONES

λ2 − a β=

λ2 +·σ 2 1 − a2

λ2 + σ 2 1 − a2

Desarrollando y simplificando se obtiene

β=

(1 + a)·λ2 λ2 + σ 2

Luego se observa que cuando t→ ∞ la estima del parámetro β no tiende al parámetro real a. ♦

♦ Solución problema 5.3 El modelo que se considera para la identificación es:

yˆ (t ) = − ay (t − 1) + bu (t − 1) En este caso:

a 

φ T (t ) = [− y (t − 1) u (t − 1)], θ =   b  

Para construir la ecuación normal

) (Φ T ·Φ )·θ = Φ T ·Y hay que hacer las siguientes multiplicaciones de matrices


64

u (1 − 1)  ·  u (i − 1)   ·  u (t − 1) tx1


SOLUCIONES

 y (1)  ·   − y (1 − 1) · − y (i − 1) · − y (t − 1)    · Φ T ·Y =  y i ( )    u (1 − 1) · u (i − 1) · u (t − 1)    ·   y (t )

tx1

que se pueden expresar en la forma: t t   2 ( 1 ) y i y (i − 1)·u (i − 1) − − ∑ ∑  i =1 · Φ T ·Φ =  t i =1 t 2  − y (i − 1)·u (i − 1) u (i − 1) ∑   ∑ i =1 i =1

  t − ∑ y (i )· y (i − 1) · Φ T ·Y =  ti =1  y (i )·u (i − 1)    ∑ i =1 Luego la ecuación normal es: t t    t  2 − − − − y ( i 1 ) y ( i 1 )· u ( i 1 ) ∑ ∑   a  − ∑ y (i )· y (i − 1) i =1  t i =1 ·  =  ti =1  t b  2 − y (i − 1)·u (i − 1)   u (i − 1) y (i )·u (i − 1)  ∑ ∑  ∑    i =1 i =1   i =1

Sustituyendo la ley de control u(t ) = − K · y (t ) en la ecuación normal se obtiene: t 1 k   a  1  2 y (i )· y (i − 1)·  − = − ( 1 ) · · y i ∑ ∑ k k 2  b  i =1 i =1    k  t

Se observa que la matriz Φ ·Φ es no singular ya que su determinante es T

1 k =0 k k2

Luego los parámetros de este modelo en lazo cerrado no se pueden identificar de forma única.

♦

65


SOLUCIONES

♦ Solución problema 5.4 a) Puesto que hay que estimar dos parámetros en el modelo la señal de entrada que se utilice en el proceso de identificación debe poseer un grado 2 de excitación persistente (EP) mayor o igual a 2. Así señales validas serían una sinusoide, una suma de sinusoides, una señal PRBS, ruido blanco, ... b) El modelo que se considera para la identificación es:

y (t ) = −θ1 · y (t − 1) − θ 2 · y (t − 2 ) + u(t − 1) En este caso:

θ1  φ (t ) = [ − y(t − 1) − y (t − 2) u(t − 1)], θ = θ 2     1  T

Para construir la ecuación normal (5.5) hay que hacer las siguientes multiplicaciones de matrices

0 u(1 − 1)   − y (1 − 1) − y (1 − 1) − y ( 2 − 1) ··· − y (t − 1)   − y ( 2 − 1) − y ( 2 − 2) u( 2 − 1)  Φ T ·Φ =  − y ( 2 − 2) ··· − y (t − 2)· 0   : : :   u(1 − 1) u( 2 − 1) ··· u(t − 1)     − y (t − 1) − y (t − 2) u(t − 1)  − y (1 − 1) − y ( 2 − 1) ··· − y (t − 1)   y (1) Φ ·Y =  0 − y ( 2 − 2 ) ··· − y (t − 2 )· :     u( 2 − 1) ··· u(t − 1)   y (t )   u(1 − 1) T

que se pueden expresar en la forma: t t t   2 − − − − y i y i y i y (i − 1)·u(i − 1)  ( 1 ) ( 1 )· ( 2 ) ∑ ∑ ∑  i =1 i =1  t i =1  t t T 2  Φ ·Φ = ∑ y (i − 1)· y (i − 2) − ∑ y (i − 2 )·u(i − 1)· y (i − 2 ) ∑  i =1  i =1 i =1 t t  t  2 − u i ( 1 ) − ∑ y (i − 1)·u(i − 1) − ∑ y (i − 2 )·u(i − 1)  ∑ i =1 i =1  i =1 

66


SOLUCIONES

  t  − ∑ y (i )· y (i − 1)    i t=1 T  Φ ·Y = − ∑ y (i )· y (i − 2 )·   i =1   t  ∑ y (i )·u(i − 1)    i =1 Luego la ecuación normal es: t t t   t   2 y i y i y i y i u i ( 1 ) ( 1 )· ( 2 ) ( 1 )· ( 1 ) − − − − − − ∑ ∑ ∑  − ∑ y(i )· y(i − 1)    i =1 i =1   θ1   i t=1  t i =1 t t 2      y(i − 1)· y(i − 2) y (i − 2 ) − ∑ y(i − 2)·u(i − 1) · θ 2 = − ∑ y(i )· y(i − 2) (1) ∑      i =1 ∑ i =1 i =1 i =1  1   t t t    t  u 2 (i − 1)  ∑ y(i )·u(i − 1)   − ∑ y(i − 1)·u(i − 1) − ∑ y(i − 2)·u(i − 1) ∑ i =1 i =1   i =1   i =1

La estima de mínimos cuadrados de los parámetros θ1 y θ2 del modelo se obtendría resolviendo el siguiente sistema de ecuaciones derivado de la ecuación normal: t t  t 2   t  y ( i − 1 ) · + y ( i − 1 )· y ( i − 2 ) · − y ( i − 1 )· u ( i − 1 ) = − y (i )· y (i − 1) θ θ ∑ ∑  1 ∑  2 ∑ i =1 i =1  i =1   i =1  t t  t   t 2  y ( i − 1 )· y ( i − 2 ) · θ + y ( i − 2 ) · θ − y ( i − 2 )· u ( i − 1 ) = − y (i )· y (i − 2 ) ∑ ∑  1 ∑  2 ∑ i =1 i =1  i =1   i =1 

(2)

Si se definen las siguientes variables t

s1 = ∑ y 2 (i − 1) i =1 t

s 5 = ∑ y 2 (i − 2 ) i =1

t

s2 = ∑ y (i − 1)· y (i − 2)

t

t

i =1

i =1

s3 = ∑ y (i − 1)·u(i − 1) s4 = ∑ y (i )· y (i − 1)

i =1 t

t

s6 = ∑ y (i − 2)·u(i − 1)

s7 = ∑ y (i )· y (i − 2)

i =1

i =1

El sistema de ecuaciones se puede expresar en la forma

s1 ·θ1 + s2 ·θ 2 − s3 = − s4 s2 ·θ1 + s5 ·θ 2 − s6 = − s7

⇒

s1 ·θ1 + s2 ·θ 2 = s3 − s4 s2 ·θ1 + s5 ·θ 2 = s6 − s7

Y resolviendo el sistema se obtienen las siguientes expresiones:

θ1 =

s5 ·( s3 − s4 ) − s2 ·( s6 − s7 ) s1 ·s5 − s22

67


θ2 =

SOLUCIONES

− s2 ·( s3 − s4 ) + s1 ·( s6 − s7 ) s1 ·s5 − s22

b) Multiplicando ambos miembros de la ecuación (2) por 1/t y haciendo la aproximación:

1 t 2 1 t 2 ·∑ x (i − 1) ≈ ·∑ x (i ) t i =1 t i =1 la ecuación normal se puede expresar en la forma: t t  t 2   t  y ( i ) · + y ( i )· y ( i − 1 ) · − y ( i )· u ( i ) = − y (i )· y (i − 1) θ θ ∑ ∑  1 ∑  2 ∑ i =1 i =1  i =1   i =1  t t  t   t 2  y ( i )· y ( i − 1 ) · θ + y ( i − 1 ) · θ − y ( i − 1 )· u ( i ) = − y (i )· y (i − 2) ∑ ∑  1 ∑  2 ∑ i =1 i =1  i =1   i =1 

Si t→ ∞ , entonces se obtiene

E[ y 2 (t )]·θ1 + E[ y (t )· y (t − 1)]·θ 2 − E[ y (t )·u(t )] = − E[ y (t )· y (t − 1)] E[ y (t )· y (t − 1)]·θ1 + E[ y 2 (t − 1)]·θ 2 − E[ y (t − 1)·u(t )] = − E[ y (t )· y (t − 2)]

(3)


E[e(t )·e(t − 1)] = 0 E[u (t )·u (t − 1)] = 0 E[e(t )·u (t )] = 0 Por otra parte, supuesto que tanto la perturbación e como la salida u son ruido blanco

E[e(i ) 2 ] = σ 2 E[u (i ) 2 ] = λ 2 Además como el valor el valor pasado de la salida y es independiente del valor actual de e o de u, se cumple:

E[ y (t − 1)·u (t )] = 0 E[ y (t − 1)·e(t )] = 0 Asimismo se consideran las siguientes aproximaciones:

68


SOLUCIONES

E[ y(t − 1) 2 ] ≈ E[ y (t − 1) 2 ] ≈ E[ y (t ) 2 ] E[u(t − 1) 2 ] ≈ E[u(t ) 2 ] = λ2 E[ e(t − 1) 2 ] ≈ E[ e(t ) 2 ] = σ 2 E[ y(t − 1)· y (t − 2)] ≈ E[ y (t )· y (t − 1)] Teniendo en cuenta toda esta información, en primer lugar, se va a calcular E[ y (t )·u (t )] :

E[ y (t )·u(t )] = E[( − a· y (t − 1) − b· y (t − 2) + u(t ) + e(t ))·u(t )] = = E [ − a· y (t − 1)·u(t )] + E[ −b· y (t − 2)·u(t )] + E[u 2 (t )] + E[ e(t )·u(t )] = = 0 + 0 + λ2 + 0 = λ2 2

En segundo lugar se va a calcular E[ y (t ) ]

E[ y (t ) 2 ] = E[( − a· y (t − 1) − b· y (t − 2) + u(t ) + e(t )) 2 ] desarrollando el cuadrado, tomando el valor esperado, y considerando únicamente los términos no nulos se obtiene la siguiente expresión

E[ y (t ) 2 ] = a 2 ·E[ y 2 (t − 1)] + 2·a ·b·E [ y (t − 1)· y (t − 2 )] + b 2 ·E[ y 2 (t − 2 )] + E[u 2 (t )] + E[ e 2 (t )] Que se puede expresar equivalentemente en la forma:

E[ y (t ) 2 ] = a 2 ·E[ y 2 (t )] + 2·a·b·E[ y (t )· y (t − 1)] + b 2 ·E[ y 2 (t )] + λ2 + σ 2 2

Despejando E [ y (t ) ] se obtiene:

E [ y (t ) 2 ] =

2·a·b·E[ y (t )· y (t − 1)] + λ2 + σ 2 1 − a2 − b2

(4)

En tercer lugar se requiere calcular E [ y (t )· y (t − 1)] :

E[ y (t )· y (t − 1)] = E[( − a· y (t − 1) − b· y (t − 2) + u(t ) + e(t ))· y (t − 1)] = = E [ − a· y 2 (t − 1)] + E[ −b · y (t − 2)· y (t − 1)] + E[u(t )· y (t − 1)] + E [ e(t )· y (t − 1)] = = − a·E[ y (t ) 2 ] − b·E[ y (t )· y (t − 1)] + 0 + 0 Despejando E [ y (t )· y (t − 1)] se obtiene:

E[ y (t )· y (t − 1)] =

69

− a·E [ y (t ) 2 ] 1+ b

(5)


SOLUCIONES

A partir de (4) y (5) es posible obtener las siguientes expresiones para E[ y (t ) ] y E [ y (t )· y (t − 1)] 2

E [ y (t ) 2 ] =

(1 + b )·(λ2 + σ 2 ) (1 − a 2 − b 2 )·(1 + b ) + 2·a 2 ·b

− a·(λ2 + σ 2 ) E[ y (t )· y (t − 1)] = (1 − a 2 − b 2 )·(1 + b) + 2·a 2 ·b En cuarto lugar se requiere calcular E [ y (t )· y (t − 2)] :

E[ y (t )· y (t − 2)] = E[( − a· y (t − 1) − b· y (t − 2) + u(t ) + e(t ))· y (t − 2)] = = E [ − a· y (t − 1)· y (t − 2)] + E [ −b · y 2 (t − 2)] + E[u(t )· y (t − 2)] + E[ e(t )· y (t − 2)] = = − a·E [ y (t ) 2 ] − b·E[ y (t )· y (t − 1)] + 0 + 0 = − a· =

(1 + b)·(λ2 + σ 2 ) − a·(λ2 + σ 2 ) · − b (1 − a 2 − b 2 )·(1 + b) + 2·a 2 ·b (1 − a 2 − b 2 )·(1 + b) + 2·a 2 ·b

− a ( λ2 + σ 2 ) − a(λ2 + σ 2 ) − a·b·(λ2 + σ 2 ) + a·b·(λ2 + σ 2 ) = E[ y (t )· y (t − 1)] = (1 − a 2 − b 2 )·(1 + b) + 2·a 2 ·b (1 − a 2 − b 2 )·(1 + b) + 2·a 2 ·b

Por lo tanto, la ecuación normal toma la forma:

E[ y 2 (t )]·θ1 + E[ y (t )· y (t − 1)]·θ 2 − E [ y (t )·u(t )] = − E[ y (t )· y (t − 1)] E[ y (t )· y (t − 1)]·θ1 + E[ y 2 (t )]·θ 2 = − E [ y (t )· y (t − 1)] Si se definen las siguientes variables

s1 = E[ y 2 (t )] =

(1 + b)·(λ2 + σ 2 ) (1 − a 2 − b 2 )·(1 + b ) + 2·a 2 ·b

s2 = E[ y (t )· y (t − 1)] =

− a·(λ2 + σ 2 ) (1 − a 2 − b 2 )·(1 + b) + 2·a 2 ·b

s3 = E[ y (t )·u(t )] = λ2 D = (1 − a 2 − b 2 )·(1 + b ) + 2·a 2 ·b El sistema de ecuaciones se puede expresar en la forma

s1 ·θ1 + s2 ·θ 2 = s3 − s2 s2 ·θ1 + s1 ·θ 2 = − s2 Y resolviendo el sistema se obtienen las siguientes expresiones:

70

(3)


θ1 =

θ2 =

SOLUCIONES

s1 ·s3 − s1 ·s2 + s22 s12 − s22

− s2 ·s3 − s1 ·s2 + s22 s12 − s22

Deshaciendo el cambio de variables y operando se obtienen las siguientes expresiones:

(1 + b )·λ2 ·[(1 − a 2 − b 2 )·(1 + b ) + 2·a 2 ·b] a + a·b + a 2 θ1 = + [(1 + b ) 2 − a 2 ](λ2 + σ 2 ) (1 + b ) 2 − a 2

θ1 =

− a·λ2 ·[(1 − a 2 − b 2 )·(1 + b ) + 2·a 2 ·b] a + a·b + a 2 + [(1 + b ) 2 − a 2 ](λ2 + σ 2 ) (1 + b ) 2 − a 2

Luego se observa que cuando t→ ∞ las estimas de los parámetros θ1 y θ2 no tienden a los parámetros reales a y b, respectivamente. ♦

71

BIBLIOGRAFIA •

ASTRÖM, K. J. y WITTENMARK, B. Sistemas Controlados por Computador. Ed. Paraninfo, 1996.

•

OGATA, K. Ingeniería de Control Moderna. Prentice Hall.1998.

•

OGATA, K. Sistemas de Control en Tiempo Discreto. Prentice Hall.1996.

72

AutomaticaII Problemas Parcial1

Recommend Documents