Lea Dfo PDF
Lea Dfo PDF
Lea Dfo PDF
PEREIRA
PROGRAMA DE TECNOLOGÍA ELÉCTRICA
Laboratorio de Investigación y Desarrollo en Electrónica
y Robótica – LÍDER
Presentado por:
iii
iv Índice general
vii
viii Lista de figuras
6.1. Simulación del filtro de Kalman usando ruido del proceso y de la medida
correlacionados. S = correlación entre el ruido del proceso y el ruido
de la medida. SF ilter = valor de S usado en el filtro de Kalman. . . . . 191
6.2. Simulación del filtro de Kalman usando ruido del proceso y de la medida
coloreados. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 200
8.1. Resultados del ejemplo que muestran una desviación estándar de los errores
de estimación de estado determinados a partir de los resultados de simulación
y determinados a partir de la matriz P del EKF. Estos resultados son para la
simulación del motor de imán permanente de dos fases. Esta tabla muestra que
la matriz P es un buen indicador de la magnitud de los errores de estimación
de estado del EKF. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 263
8.2. Ejemplo [] estimación de errores para el filtro de Kalman extendido,
el filtro de Kalman unscented estándar con 2n puntos sigma y el filtro
de Kalman unscented esférico con (n + 1) puntos sigma. Generalmen-
te el UKF estándar se desempeña mejor. El desempeño y el esfuerzo
computacional del UKF esférico se encuentran entre los de el EKF y
el UKF estándar . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 300
ix
Lista de algoritmos
xi
xii
Capı́tulo 1
1.1. Introducción
E n este capı́tulo se revisan algunos aspectos básicos del álgebra lineal. La teorı́a de estima-
ción óptima de estados depende fuertemente de la teorı́a de matrices, incluyendo cálculo
matricial; de manera que la primera parte está dedicada a una revisión de este tópico. La
estimación óptima de estados se aplica tanto a sistemas lineales como a no lineales, aunque
es más directa para los lineales. Por lo tanto, se dará una breve revisión de estos sistemas.
Las últimas secciones se dedicarán a revisar los conceptos de estabilidad, controlabilidad y
observabilidad de los sistemas lineales. Estos conceptos son necesarios para comprender el
material de estimación óptima de estados que se estudiará más adelante.
[1 5 − 3 0]
es un vector de 4 elementos. Este vector se denomina 1 × 4 porque tiene una fila y cuatro
columnas, también se denomina vector fila. El vector
[ ]⊤
x = cos ωk sen ωk − sen ωk cos ωk 0 1
1
2 Capı́tulo 1. Elementos de cálculo matricial
es una matriz 3 × 4 puesto que tiene 3 filas y cuatro columnas. El número de filas y columnas
se define como la dimensión de la matriz. Se define el rango de una matriz como el número
de filas o columnas linealmente independientes [Kreyszig, 2011]. El rango de la matriz A se
indica como ρ(A). El rango de una matriz siempre es menor o igual que el número de filas
o columnas, según el caso que se contemple. Por ejemplo, la matriz
[ ]
2 3 1 −2
A=
4 6 2 −4
es uno (Matlab: rank(A)) debido a que tiene solamente una fila linealmente independiente;
las dos filas son múltiples entre sı́. Por otra parte, la matriz
1 2 5
A = 3 4 −1
−1 3 5
tiene rango 3 puesto que todas sus filas y columnas son linealmente independientes. Esto es,
existen escalares α1 , α2 y α3 tales que
[ ] [ ] [ ] [ ]
α1 1 2 5 + α2 3 4 −1 + α3 −1 3 5 = 0 0 0
de modo que las tres filas son linealmente independientes. Lo mismo para las columnas. Esto
es, hay escalares α1 , α2 y α3 diferentes de cero tales que
1 2 5 0
α1 3 + α2 4 + α3 −1 = 0
−1 3 5 0
de modo que las tres columnas son linealmente independientes. Una matriz cuyos elementos
son todos cero tiene un rango de cero. Una matriz n × m cuyo rango es igual al min(n, m)
es de rango completo. La nulidad de una matriz A n × m es igual a [m − ρ(A)].
La transpuesta de una matriz o vector se puede obtener por el cambio de todas las filas
a columnas, y todas las columnas a filas o de columna a fila en el caso de un vector. Por
1.3. Álgebra de matrices 3
entonces A⊤ es la matriz p × r
a11 · · · ar1
..
A⊤ = ... ..
. .
a1p · · · arp
entonces AH es la matriz
2 1 + j2 −j3
AH = 3 + j4 2 5 + j2
−j 1 − j3 5
donde aji denota la conjugada compleja. Como resultado de esta definición, los elementos
de la diagonal aii de una matriz hermı́tica son números reales (ya que aii = aii ), mientras
que otros elementos pueden ser complejos.
Como un ejemplo de matrices hermı́ticas 2 × 2 se tienen las matrices de Pauli [Nakahara,
2003]:
[ ] [ ] [ ]
0 1 0 −i 1 0
σ1 = σ2 = σ3 =
1 0 i 0 0 −1
C = AB := (cij )m×p
es decir:
a11 ··· a1n b11 ··· b1p
.. .. .. .. .. ..
C = Am×n Bn×p = . . . × . . .
am1 ··· amn m×n
bn1 ··· bnp n×p
a11 b11 + · · · + a1n bn1 ··· a11 b1p + · · · + a1n bnp
.. .. ..
= . . .
am1 b11 + · · · + amn bn1 ··· am1 b1p + · · · + amn bnp m×p
sigue:
x1
[ ]
x2
x⊤ x = x1 x2 · · · xn . = x21 + x22 + · · · + x2n (1.4)
..
xn
x1 x21 x1 x2 ··· x1 xn
x2 [
] x2 x1 x22 ··· x2 xn
xx⊤ = . x1 x2 ··· xn = . .. .. .. (1.5)
.. .. . . .
xn xn x1 xn x2 ··· x2n
Supóngase que se tiene una matriz H(p×n) y una matriz P(n×n) . Entonces H⊤ es una
matriz n × p, y se puede calcular el producto matricial HPH⊤
(p×p) [Poznyak, 2008].
h11 · · · h1n p11 · · · p1n h11 ··· hp1
⊤ .. .. .. .. . .
.. ... .. ..
HPH = . . . . .. . .
hp1 · · · hpn pn1 · · · pnn h1n ··· hpn
∑ ∑
h1j pjk h1k · · · h1j pjk hpk
j,k j,k
.. .. ..
= . . . (1.6)
∑ ∑
hpj pjk h1k · · · hpj pjk hpk
j,k j,k
En general:
AB ̸= BA (1.8)
2. Asociativa:
5. Inversa:
En general
6. Transpuesta:
(A + B)⊤ = A⊤ + B⊤ (1.14)
(AB)⊤ = B⊤ A⊤ (1.15)
En general
7. Conjugada Hermı́tica:
1.3. Álgebra de matrices 7
(A + B)H = AH + BH (1.18)
(AB)H = BH AH (1.19)
8. Traza:
∑ ∑
tr(A) = aii = λi ; donde λi = eig(A) (1.22)
i i
9. Determinante:
Sea A una matriz n × n. Usando la notación Aij para representar la matriz formada
borrando la i-ésima fila y la j -ésima columna de A. Entonces, el determinante de A
se define como
∑
n
|A| = (−1)i+j aij |Aij | (1.26)
j=1
para cualquier valor de j ∈ [1, n]. Ambos desarrollos conducen a un mismo resultado.
A partir de la definición se puede observar para los casos (2 × 2) y (3 × 3):
[ ]
a a12
det(A) = det 11 = a11 a22 − a12 a21 (1.28)
a21 a22
a11 a12 a13
det(A) = det a21 a22 a23 =
a31 a32 a33
[ ] [ ] [ ]
a a23 a a23 a a22
= a11 det 22 − a12 det 21 + a13 det 21 =
a32 a33 a31 a33 a31 a32
= a11 (a22 a33 − a23 a32 ) − a12 (a21 a33 − a23 a31 ) + a13 (a21 a32 − a22 a31 )
(1.29)
Propiedades:
Las matrices cuadradas se pueden multiplicar por sı́ mismas en forma repetida del
mismo modo que los números ordinarios, puesto que ellas tienen siempre el mismo
número de filas y columnas. Esta multiplicación repetida se puede describir como la
“potencia de la matriz”, un caso especial del producto de matrices. Las matrices “rec-
tangulares” no tienen el mismo número de filas y columnas por lo que no se pueden
elevar a una potencia. Una matriz An×n elevada a un entero positivo p se define como
· · · A}
Ap = |AA{z
p veces
A0 = I
I es la matriz identidad.
(ii) Multiplicación por un escalar:
(αA)p = αp Ap
Norma euclidiana
La norma dos de un vector columna de números reales, también llamada la norma euclidiana,
se define como sigue:
√ √
||x||2 = x⊤ x = x21 + · · · + x2n (1.36)
Esto significa que la norma euclidiana al cuadrado equivale a la traza del producto del vector
por su transpuesta, es decir:
( )
x⊤ x = Tr xx⊤ (1.38)
Dado que las matrices cuadradas 2 × 2 son de uso muy frecuente en la práctica, es de interés
analizar este caso ya que posee algunas propiedades de interés, las cuales permiten simplificar
los cálculos. Considérese la matriz
[ ]
a11 a12
A= (1.39)
a21 a22
10 Capı́tulo 1. Elementos de cálculo matricial
1. Determinante y traza.
2. Valores propios
[ ] [ ] [ ]
λ 0 a a12 λ − a11 −a12
|[λI − A]| = 0 = − 11 =
0 λ a21 a22 −a21 λ − a22
(λ − a11 )(λ − a22 ) − a12 a21 = λ2 − (a11 + a22 )λ + a11 a22 − a12 a21
Entonces
λ2 − λ · tr(A) + det(A) = 0, (1.42)
1 [ √ ] 1[ √ ]
λ1 = tr(A) + tr(A)2 − 4det(A) , λ2 = tr(A) − tr(A)2 − 4det(A) , (1.43)
2 2
λ1 + λ2 = tr(A), λ1 λ2 = det(A). (1.44)
3. Autovectores
[ ]
a12
v1 ∝ (1.45)
λ1 − a11
[ ]
a12
v2 ∝ (1.46)
λ2 − a11
4. Inversa
[ ]
−1 1 a22 −a12
A = (1.47)
det(A) −a21 a11
Ejemplo 1.1
Realimentación serie–serie con un elemento activo (NFC) [Avendaño, 2007].
El circuito de la Fig. 1.1(a), constituye una red lineal con realimentación a tra-
vés de un circuito activo. Se puede definir la red conformada por Ra , Rb y el
amplificador Ao1, como una fuente de tensión controlada por tensión, es decir,
( )
Rb .
vβ = 1 + vm = µvm (1.48)
Ra
La red equivalente se muestra en la Fig. 1.1(b), donde se indica la relación definida
en la ecuación (1.48). Las resistencias de entrada y salida de Ao1 y Ao2 están
dadas por Ria , Roa , Rib y Rob , respectivamente.
1.3. Álgebra de matrices 11
Entonces:
det(Z) = zi zo − βz µz
Tr(Z) = zi + zo ≈ Ria
[ √ ]
1[ √ ] R
ia Rc
λ1 = zi + zo + (zi − zo ) + 4βz µz ≈
2 1 + 1 − µAo
2 2 Ria
[ √ ]
1[ √ ] R
ia R c
λ2 = zi + zo − (zi − zo )2 + 4βz µz ≈ 1 − 1 − µAo
2 2 Ria
Ejemplo 1.2
Realimentación Serie–Paralelo [Avendaño, 2007].
El modelo de realimentación serie–paralelo se muestra esquemáticamente en la
Figura 1.2.
Es decir,
[ ] [ ]
vs i
=H s
0 vo
donde
[ ]
zi β
H= (1.51)
µ yo
donde
[ ( )]
µ
|H| = yo zi − βµ = yo zi 1+β −
yo zi
o
|H| = yo zi (1 + βAva )
con
µ
Ava = − .
yo z i
De esta expresión se pueden calcular los parámetros de transferencia, v.gr.:
Ganancia de tensión,
vo µ −µ Ava
Av = =− =− =
vs |H| yo zi (1 + βAva ) 1 + βAva
En este problema, puesto que el modelo de la matriz tiene dimensiones hı́bridas,
es posible calcular algunas funciones de la misma, tales como el determinante y
la matriz inversa; sin embargo, la traza no está definida ya que se tendrı́a la suma
de [Ω] + [S] (ver Ecuación (1.51)), lo cual no es posible. Esto se debe tener en
cuenta en los sistemas fı́sicos, pues ocurre con alguna frecuencia.
En ella A y BD−1 C son matrices de dimensión n×n, mientras que D es una matriz cuadrada
que puede tener una dimensión r × r diferente, en consecuencia B tendrá dimensión n × r y
C dimensión r × n. Los nombres alternativos para esta fórmula son el lema de inversión de
matrices, la fórmula de Sherman-Morrison-Woodbury [Sherman and Morrison, 1950] o sólo
la fórmula de Woodbury [Woodbury, 1950]. Sin embargo, la identidad apareció en varios
artı́culos antes del informe de Woodbury.
Demostración:
Se construye una matriz aumentada A, B, C y D y su correspondiente inversa:
[ ]−1 [ ]
A B E F
= (1.53)
C D G H
Las submatrices en (1.54) y (1.55) se abren para formar ocho ecuaciones matriciales:
AE + BG = I (1.56)
AF + BH = 0 (1.57)
CE + DG = 0 (1.58)
CF + DH = I (1.59)
EA + FC = I (1.60)
EB + FD = 0 (1.61)
GA + HC = 0 (1.62)
GB + HD = I (1.63)
De (1.56) y (1.58)
AE = I − BG
DG = −CE
G = −D−1 CE
AE = I + BD−1 CE
AE − BD−1 CE = I
(A − BD−1 C)E = I
E = (A − BD−1 C)−1 (1.64)
1.3. Álgebra de matrices 15
De (1.61) y (1.64)
EB + FD = 0
FD = −EB
F = −(A − BD−1 C)−1 BD−1 (1.65)
De (1.58) y (1.64)
DG = −CE
G = −D−1 C(A − BD−1 C)−1 (1.66)
De (1.59) y (1.65)
De (1.63) y (1.62)
HD = I − GB
GA = −HC
G = −HCA−1
HD = I + HCA−1 B
H(D − CA−1 B) = I
H = (D − CA−1 B)−1 (1.68)
De (1.62) y (1.68)
De (1.57) y (1.68)
AF = −BH
AF = −B(D − CA−1 B)−1
F = −A−1 B(D − CA−1 B)−1 (1.70)
De (1.60) y (1.70)
EA = I − FC
EA = I + A−1 B(D − CA−1 B)−1 C
E = A−1 + A−1 B(D − CA−1 B)−1 CA−1 (1.71)
La prueba se completa mediante la combinación de, v.gr., (1.64) y (1.71), con la condición
16 Capı́tulo 1. Elementos de cálculo matricial
Aq = λq (1.73)
Por lo tanto
(A − λI)q = 0 (1.74)
Para el análisis de valores propios se procede como sigue para satisfacer esta ecuación o bien
q = 0, lo cual es interesante, o la matriz A − λI debe reducir q al vector nulo (un único
punto). Para que esto suceda A − λI debe ser singular. Por lo tanto,
Por lo tanto, el análisis de valores propios procede ası́ (i ) se resuelve la ecuación anterior
para encontrar los autovalores λi y entonces (ii) se sustituyen en la Ecuación (1.73) para
encontrar los autovectores. Por ejemplo, si
[ ]
4 −5
A= (1.76)
2 −3
entonces
λ2 − λ − 2 = 0
(1.78)
(λ + 1)(λ − 2) = 0
1.3.7. Diagonalización
Si se colocan los autovectores en las columnas de una matriz
[ ]
Q = q1 q2 · · · qd (1.82)
AQ = QΛ (1.84)
Q−1 AQ = Λ (1.85)
esto muestra que cualquier matriz cuadrada se puede convertir en una forma diagonal (siem-
pre que tenga valores propios distintos; ver [Strang, 2006]). A veces no habrá valores propios
d distintos y en ocasiones ellos serán complejos.
Q⊤ AQ = Λ (1.86)
A = QΛQ⊤ (1.87)
18 Capı́tulo 1. Elementos de cálculo matricial
lo cual se conoce como el teorema espectral. Éste dice que cualquier matriz real simétrica
puede representarse como la anterior donde las columnas de Q contienen los autovectores y
Λ es una matriz diagonal que contiene los autovalores, Λi . Equivalentemente
λ1 q1
[
] λ
2 q2
A = q1 q2 · · · qd . . (1.88)
.. ..
λd qd
∑
d
A= λk qk q⊤
k (1.89)
k=1
1. Conmutativo
A ◦ B = B ◦ A, (1.91)
2. Lineal
3. Asociativo
A ◦ (B ◦ C) = (A ◦ B) ◦ C, (1.93)
1.3. Álgebra de matrices 19
4. Distributivo
A ◦ (B + C) = A ◦ B + A ◦ C. (1.94)
5. Rango.
6. Matrices diagonales.
Ejemplo 1.3
Verificar las relaciones de los productos mixtos de Hadamard de la Ecuación (1.96),
para los valores de A, B, D, E como se especifican a continuación:
[ ] [ ] [ ] [ ]
1 0 0
1 2 −1 −2 4 1 −2 0
A= , B= , D= , E = 0 −1 0 .
−1 3 −2 1 −3 −1 0 1
0 0 4
Dadas dos matrices, A(m×n) y B(p×q) , el producto de Kronecker, producto directo o pro-
ducto tensorial se define como la matriz A ⊗ B(mp×nq) (Matlab: kron(A, B)) [Pedersen and
20 Capı́tulo 1. Elementos de cálculo matricial
Petersen, 2012].
a11 B a12 B ··· a1n B
a21 B a22 B ··· a2n B
A⊗B= . .. .. .. (1.97)
.. . . .
am1 B am2 B · · · amn B
A ⊗ B ̸= B ⊗ A, en general. (1.98)
2. Asociatividad.
(A ⊗ B) ⊗ C = A ⊗ (B ⊗ C). (1.99)
4. Transpuesta
(A ⊗ B)⊤ = A⊤ ⊗ B⊤ . (1.101)
5. Distributividad
A ⊗ (B + C) = A ⊗ B + A ⊗ C. (1.102)
6. Rango
7. Escala
8. Traza
9. Inversa.
1.3. Álgebra de matrices 21
10. Determinante
12. Ortogonalidad.
Nótese que el producto de Kronecker está definido independientemente del orden de las ma-
trices involucradas, por lo que es un concepto más general que la multiplicación de matrices.
A continuación se hará la demostración de algunas de las propiedades relacionadas anterior-
mente. Dichas demostraciones tienen el carácter de teoremas y en esa forma se plantean.
Teorema 1.1
Sea A ∈ Rm×n y B ∈ Rp×q , entonces
(A ⊗ B)−1 = A−1 ⊗ B−1 .
Demostración: De la Ecuación (1.100):
(A ⊗ B)(C ⊗ D) = AC ⊗ BD.
Sustituyendo, se obtiene sucesivamente:
(A ⊗ B)(A ⊗ B)−1 = (A ⊗ B)(A−1 ⊗ B−1 )
= AA−1 ⊗ BB−1
= Im ⊗ In = Imn
22 Capı́tulo 1. Elementos de cálculo matricial
Teorema 1.2
Sea A ∈ Rm×n y B ∈ Rp×q , entonces
A ⊗ B = (A ⊗ Ip )(In ⊗ B) = (Im ⊗ B)(A ⊗ Iq ). (1.110)
Demostración: De acuerdo a la definición del producto de Kronecker y a la
multiplicación de matrices, se tiene
a11 B a12 B · · · a1n B
a21 B a22 B · · · a2n B
A ⊗ B = ... .. .. =
.
..
. .
am1 B am2 B · · · amn B
a11 Ip a12 Ip · · · a1n Ip B 0 · · · 0
a21 Ip a22 Ip · · · a2n Ip 0 B · · · 0
= ... .. .. . . . .
. .. .. . . ..
..
. .
am1 Ip am2 Ip · · · amn Ip 0 0 · · · B
= (A ⊗ Ip )(In ⊗ B),
a11 B a12 B · · · a1n B
a21 B a22 B · · · a2n B
A ⊗ B =
... .. .. =
.
..
. .
am1 B am2 B · · · amn B
B 0 · · · 0 a11 Iq a12 Iq · · · a1n Iq
0 B · · · 0 a21 Iq a22 Iq · · · a2n Iq
=
... ... . . .. ..
. .. .. .
.. ..
. .
0 0 · · · B am1 Iq am2 Iq · · · amn Iq
A ⊗ B = (Im ⊗ B)(A ⊗ Iq ).
Corolario 1.1
Sea A ∈ Rm×m y B ∈ Rn×n , entonces
A ⊗ B = (A ⊗ In )(Im ⊗ B) = (Im ⊗ B)(A ⊗ In ). (1.111)
Esto significa que Im ⊗ B y A ⊗ In son conmutativos para matrices cuadradas
tales como A y B.
1.3. Álgebra de matrices 23
A partir del Teorema 1.2, se puede demostrar el siguiente teorema de productos mixtos.
Teorema 1.3
Sea A ∈ Rm×n , B ∈ Rq×r , C ∈ Rn×p y D ∈ Rr×s , entonces
(A ⊗ B)(C ⊗ D) = (AC) ⊗ (BD). (1.112)
Demostración: De acuerdo al Teorema 1.2, se tiene
(A ⊗ B)(C ⊗ D)
= (A ⊗ Iq )(In ⊗ B)(C ⊗ Ir )(Ip ⊗ D)
= (A ⊗ Iq )[(In ⊗ B)(C ⊗ Ir )](Ip ⊗ D)
= (A ⊗ Iq )(C ⊗ B)(Ip ⊗ D)
= (A ⊗ Iq )[(C ⊗ Iq )(Ip ⊗ B)](Ip ⊗ D)
= [(A ⊗ Iq )(C ⊗ Iq )] [(Ip ⊗ B)(Ip ⊗ D)]
= [(AC ⊗ Iq )] [(Ip ⊗ BD)]
= (AC) ⊗ (BD).
Ejemplo 1.4
Sean las matrices
[ ] [ ] [ ] [ ]
1 2
1 2 3 2 3 1 0 1 0
A= , B= , C = −1 1 , D=
4 5 6 1 2 1 1 −1 2
−2 1
donde m = 2, n = 3, k = 2, l = 2, q = 2, s = 4. Verificar las relaciones de la
Ecuación (1.108).
Solución:
La relación directa está dada por
1 0 1 0 2 0 2 0 3 0 3 0
1 1 −1 2 2 2 −2 4 3 3 −3 6
A⊗D=
4 0 4 0 5 0 5 0 6 0 6 0
4 4 −4 8 5 5 −5 10 6 6 −6 12
1 0 1 0 0 0 0 0 0 0 0 0
1 1 −1 2 0 0 0 0 0 0 0 0
0 0 0 0 1 0 1 0 0 0 0 0
I3 ⊗ D =
0 0 0 0 1 1 −1 2 0 0 0 0
0 0 0 0 0 0 0 0 1 0 1 0
0 0 0 0 0 0 0 0 1 1 −1 2
1. Producto
2. Traza
3. Suma
4. Escala
5. Producto matriz–vector
A continuación se demostrará la relación que existe entre matrices y el operador vec denomi-
nada vectorización, la cual permite la transformación de un producto matricial en un vector.
Este resultado será utilizados más adelante en el cálculo de las ecuaciones de Sylvester y
Lyapunov.
Teorema 1.4
Sean las matrices A ∈ Rm×n , B ∈ Rp×q y X ∈ Rn×p . Demostrar que se cumple
la relación (Ecuación (1.113)):
Entonces,
( )
(AXB).k = b⊤
.k ⊗ A vec(X)
26 Capı́tulo 1. Elementos de cálculo matricial
es decir,
Ejemplo 1.5
Verificar la relación (1.113) a partir de las siguientes matrices:
[ ] [ ] [ ]
1 2
1 2 3 2 3 4 1
A= , B= , X= 1 1 .
4 3 1 1 3 2 1
2×3 2×4 −1 1 3×2
Solución:
1
2
[ ] 12
1 12 2 4 6
vec(AXB) = vec =
2
2 6 4 2
4
4
2
2 4 6 1 2 3 1
−2 6 2 −1 3 1 1 2
3 3 6 9 1 12
6 9
−3 9 3 −3 9 3 −1 6
(B⊤ ⊗ A)vec(X) =
4 2 2
=
8 12 2 4 6
−4 12 4 −2 6 2 1 4
1 2 3 1 2 3 1 4
−1 3 1 −1 3 1 2
Con lo cual queda demostrada la relación pedida.
1.3. Álgebra de matrices 27
Ejemplo 1.6
A continuación se definen algunas cantidades vectoriales y matriciales. Verificar
la Ecuación (1.117).
[ ] [ ] [ ] [ ]
1 −2 1 −2
2 3
a= 2 , c= 1 , B= , X = −1 3
1 −4
3 5 2 −1
[ ][ ] [ ] [−2]
1 −2
2 3 1 −1 2
a⊤ XBX⊤ c = [1 2 3] −1 3 1 = 99
1 −4 −2 3 −1
2 −1 5
Nota sobre la función vec y Matlab. Matlab no tiene la función vec, pero posee
una función llamada reshape, la cual puede usarse para este propósito. Especı́ficamente, si
A ∈ Rm×n , entonces se puede obtener el vector a con el siguiente enunciado de Matlab:
a = reshape(A, m*n, 1);
Existe una forma más corta para realizar la misma tarea, esto es,
a = A(:);
la cual es similar a la expresión anterior. Para reconstruir la forma de matriz, es decir, de
a → A, se puede usar de nuevo la función reshape:
A = reshape(a, p, q); donde pq = mn, es decir, mn|p y mn|q (Resto = 0).
Ejemplo 1.7
Sea A ∈ R3×2 , dada por:
1
[ ] 3
1 2
1
A = 3 −1 ⇒ [m, n] = size(A) = [3, 2]; a = reshape(A, m ∗ n, 1) =
2
1 4 −1
4
[ ]
1 1 −1
A1 = reshape(a, n, m) ⇒ A1 = .
3 2 4
A ⊕ B = A ⊗ In + Im ⊗ B. (1.118)
⊕
N
Ai = A1 ⊕ A2 ⊕ · · · ⊕ AN . (1.120)
i=0
2. Rango
( )
⊕
N ∑
N
ρ Ai = ρ(Ai ). (1.122)
i=0 i=0
3. Determinante
Si cada Ai es una matriz cuadrada, entonces
( N )
⊕ ∏N
det Ai = det(Ai ). (1.123)
i=0 i=0
1.3. Álgebra de matrices 29
4. Ortogonalidad
Si A1 y A2 son ortogonales, entonces A1 ⊕ A2 es ortogonal.
5. Exponencial
Sea A ∈ Rm×m y B ∈ Rn×n , entonces
Los tres resultados son idénticos como era de esperarse. Recuérdese que en Matlab
la exponencial de una matriz se define como expm(A).
Teorema 1.5
Teorema de estabilidad de Lyapunov.
Sea V una función no negativa sobre Rm y sea V̇ la derivada en el tiempo de V
a lo largo de las trayectorias de la dinámica del sistema (1.126):
dV dV dx dV
V̇ = = = ẋ (1.127)
dt dx dt dx
Sea Br = Br (0) una bola de radio r alrededor del origen. Si ∃ r > 0 tal que V sea
definida positiva y V̇ sea semidefinida negativa para todo x ∈ Br , entonces x = 0
es localmente estable en el sentido de Lyapunov. Si V es definida positiva y V̇ es
definida negativa en Br , entonces x = 0 es localmente asintóticamente estable.
Si V satisface una de las condiciones anteriores, se dice que V es una función (local) de
Lyapunov para el sistema. Siempre es posible encontrar una función de Lyapunov para un
sistema lineal de la forma descrita por (1.126), sin perturbación externa, esto es, con u = 0.
En este caso, se puede escoger una función candidata de Lyapunov de la forma cuadrática:
V = x⊤ Px (1.128)
V̇ = ẋ⊤ Px + x⊤ Pẋ
= (Ax)⊤ Px + x⊤ P(Ax)
= x⊤ A⊤ Px + x⊤ PAx
= x⊤ (A⊤ P + PA)x
= −x⊤ Qx (1.129)
El sistema matricial de ecuaciones (1.129), tiene la forma general de un sistema del tipo
AX + XB = C. Este sistema general debe ser adecuado para el sistema particular de Lya-
punov, haciendo B = A⊤ y C = −Q. Si la matriz Q es definida positiva, entonces el sistema
es asintóticamente estable.
A continuación, se abordará la solución del sistema matricial planteado por la aplicación
del método de Lyapunov, dado por
A⊤ P + PA = −Q. (1.130)
AX + XB = C, (1.131)
Despejando vec(X):
AX + XA = C (1.137)
⊤ −1
vec(X) = (I ⊗ A + A ⊗ I) vec(C); ⇒ vec(X) → X ∈ R m×n
. (1.138)
% AX + XB = C
x1=(kron(I3,A)+kron(B’,I2))\C(:);
X1=vec2mat(x1,2)’;
% Generalized Sylvester
% AXD + EXB = C
x2=(kron(D’,A)+kron(B’,E))\C(:);
X2=vec2mat(x2,2)’;
% Lyapunov
% AX + XA = C1
x3=(kron(I2,A)+kron(A’,I2))\C1(:);
X3=vec2mat(x3,2)’;
Resultados:
−1.6766
1.4203 [ ]
−0.1206 −1.6766 −0.1206 0.8241
x1 = =⇒ X1 =
−0.1139 1.4203 −0.1139 −0.1413
0.8241
−0.1413
−0.3531
0.3959 [ ]
−0.0460 −0.3531 −0.0460 0.1516
x2 = =⇒ X2 =
−0.0283 0.3959 −0.0283 0.0288
0.1516
0.0288
4.50 [ ]
−1.10 4.50 −2.35
x3 =
−2.35
=⇒ X3 =
−1.10 0.75
0.75
solución del sistema perturbado x̂ de la solución del sistema x. O sea, saber de antemano
cual es la sensibilidad de la solución a las perturbaciones en los datos de entrada. Para tales
fines existe una cantidad que da una medida de esa sensibilidad. Dicha cantidad se denomina
número de condición [Householder, 1975]. El número de condición se define por
||Ax||
||A|| = sup = sup ||Ax|| (1.140)
x∈Rn ||x|| ||x||=1
Se hará el análisis por partes, primero suponiendo que existe un error de redondeo en b y
que A es exacta, luego recı́procamente, para finalmente hacer el caso conjunto.
Sea la ecuación lineal dada por
Ax = b (1.141)
A(x + ∆x) = b + ∆b
Ax + A∆x = b + ∆b
A∆x = ∆b
∆x = A−1 ∆b
||∆x|| = ||A−1 · ∆b|| ≤ ||A−1 || · ||∆b||
y
1 1
≤ ||A|| ·
||x|| ||b||
Donde κ(A) = ||A|| · ||A−1 ||, es el número de condición, como se habı́a definido antes.
1.4. Número de condición 35
O sea
||∆x|| ||∆b||
≤ κ(A) (1.142)
||x|| ||b||
(A + ∆A)(x + ∆x) = b
∆A · (x + ∆x) + A · ∆x = 0
∆A · x + (A + ∆A)∆x = 0
∆x = −(A + ∆A)−1 · ∆A · x (1.143)
Entonces
||A−1 ||
||(A + ∆A)−1 || ≤ (1.144)
1 − ||A−1 || · ||∆A||
Es decir,
||∆x||
≤ ||(A + ∆A)−1 ∆A||
||x||
≤ ||(A + ∆A)−1 || · ||∆A|| (1.146)
o sea,
||∆A||
κ(A) ·
||∆x|| ||A||
= (1.147)
||x|| ||∆A||
1 − κ(A) ·
||A||
(A + ∆A)(x + ∆x) = b + ∆b
(A + ∆A)∆x = ∆b − ∆A · x
∆x = (A + ∆A)−1 (∆b − ∆A · x)
Por lo tanto,
( )
||∆x|| κ(A) ||∆b|| ||∆A||
≤ · + (1.148)
||x|| ||∆A|| ||b|| ||A||
1 − κ(A) ·
||A||
Ejemplo 1.10
Sea el siguiente sistema de la forma Ax = b, donde
10 7 8 7 32 0.1
7 5 6 5 23 −0.1
A=
6 10 9
b = , con ∆b =
0.1
, .
8 33
7 5 9 10 31 −0.1
1.5. Cálculo vectorial y matricial 37
Supóngase que la matriz A tiene elementos que son funciones del tiempo. Se sabe que
AA−1 = I; esto es, AA−1 es una matriz constante y por lo tanto tiene derivada en el
tiempo igual a cero, entonces la derivada en el tiempo de AA−1 se puede calcular como
d d
(AA−1 ) = ȦA−1 + A (A−1 ) (1.151)
dt dt
d(A−1 )
Puesto que esto es cero, puede resolverse para dt como
d
(A−1 ) = −A−1 ȦA−1 (1.152)
dt
Nótese que para el caso especial de una función escalar A(t), esto se reduce a la ecuación
familiar
( )
d 1 ∂(1/A(t)) dA(t)
=
dt A(t) ∂A(t) dt
Ȧ(t)
=− (1.153)
A2 (t)
Con estas definiciones se puede calcular la derivada parcial del producto escalar de dos
vectores. Sean x e y vectores columna {x, y ∈ Rn }. Entonces
x⊤ y = x1 y1 + · · · + xn yn (1.154)
∂(x⊤ Ax)
= 2x⊤ A (1.159)
∂x
Entonces
∂f1 (x) ∂f1 (x)
···
∂x1 ∂xn
∂f (x)
=
..
.
..
.
..
.
(1.160)
∂x ∂f (x)
m ∂f (x)
m
···
∂x1 ∂xn
Bien sea que f (x) o x sea transpuesta, entonces la derivada parcial también es transpuesta.
( )⊤
∂f ⊤ (x) ∂f (x)
=
∂(x) ∂(x)
( )⊤
∂f (x) ∂f (x) (1.161)
=
∂(x)⊤ ∂(x)
∂f ⊤ (x) ∂f (x)
=
∂(x)⊤ ∂(x)
Con estas relaciones, se pueden derivar las siguientes igualdades. Supóngase que A ∈ Rm×n
y x ∈ Rn×1 . Entonces
∂(Ax)
=A
∂x (1.162)
∂(x⊤ A⊤ )
= A⊤
∂x
Por lo tanto
∂Tr(ABA⊤ )
= AB⊤ + AB (1.165)
∂A
Si B es simétrica entonces
∂Tr(ABA⊤ )
= 2AB si se cumple que B = B⊤ (1.166)
∂A
1.6. Resumen
En este capı́tulo se ha hecho una revisión de aspectos básicos de álgebra lineal incluyendo
una introducción a productos de Hadamard y de Kronecker, los cuales no son muy corrientes
en los cursos introductorios. La razón es que se presentan algunas procesos en la práctica que
requieren esas aplicaciones, tales como en la solución de la ecuación de Sylvester o Lyapunov
que tiene ocurrencia en sistemas de control. Esto conduce a un conjunto de resultados que se
han expuesto en el texto en forma de teoremas con el objeto de darles formalidad y también
como ejercicio para el lector. También se ha hecho un estudio del concepto de número de
condición. Esto es de utilidad en el proceso de análisis de estabilidad numérica de algunos
sistemas dinámicos. En la segunda parte se abordó el problema del cálculo de matrices,
procedimiento que conduce a la obtención de resultados directos para diversas aplicaciones
de matrices y vectores. Se ha procurado realizar los procedimientos por el método de paso
a paso con el fin de que el lector pueda seguir fácilmente las demostraciones. Finalmente, se
hace una pequeña introducción a los sistemas lineales y no lineales y a la discretización de
los procesos dinámicos.
Bibliografı́a
43
Capı́tulo 2
Elementos de sistemas
dinámicos
2.1. Introducción
45
46 Capı́tulo 2. Elementos de sistemas dinámicos
o lo que es lo mismo,
∫ t
−At −At0
e x(t) − e x(t0 ) = e−Aτ Bu(τ )dτ.
t0
2.3. Sistemas no lineales 47
o sea,
∫ t
−At −At0
e x(t) = e x(t0 ) + e−Aτ Bu(τ )dτ.
t0
por esta razón, eAt se denomina matriz de transición de estado del sistema.
Es la matriz que describe como cambia el estado de su condición inicial en la
ausencia de entradas externas. (Matlab: expm(A) ̸= exp(A), exp(A) calcula la
exponencial de una matriz, elemento por elemento).
ẋ = f (x, u, v)
(2.5)
y = h(x, w)
donde f (·) y h(·) son funciones vectoriales arbitrarias [Khalil, 1996]. Se utiliza
v para indicar el ruido del proceso y w para indicar el ruido de la medición. Si
f (·) y h(·) son funciones explı́citas del tiempo, entonces el sistema es variante en
el tiempo. De otra forma, el sistema es invariante en el tiempo. Si f (x, u, v) =
Ax + Bu + v y h(x, w) = Hx + w, entonces el sistema es lineal. De otra
forma es no lineal. Con el fin de aplicar herramientas a partir de la teorı́a
de los sistemas lineales a los sistemas no lineales, se necesita linealizar estos
sistemas. En otras palabras, se necesita encontrar un sistema lineal que sea
aproximadamente igual al sistema no lineal. Para ver como se hace esto, se
empieza con una función no lineal f (x) de un escalar x. Se expande f (x) en
una serie Taylor alrededor de un punto de operación nominal (también llamado
punto de linealización) x = x̄. Definiendo x̃ = x − x̄:
∂f 1 ∂ 2 f 2 1 ∂ 3 f 3
f (x) = f (x̄) + x̃ + x̃ + x̃ + · · · (2.6)
∂x x̄ 2! ∂x2 x̄ 3! ∂x3 x̄
Taylor como
( ) ( )2
∂ ∂ 1 ∂ ∂
f (x) = f (x̄) + x̃1 + · · · + x̃n f (x) + x̃1 + · · · + x̃n f (x)
∂x1 ∂xn x̄ 2! ∂x1 ∂xn
( )3 x̄
1 ∂ ∂
+ x̃1 + · · · + x̃n f (x) + · · · (2.7)
3! ∂x1 ∂xn
x̄
operación nominal (x̄, ū, v̄), obteniéndose una aproximación lineal como sigue;
ẋ = f (x, u, v)
∂f ∂f ∂f
≈ f (x̄, ū, v̄) + (x − x̄) + (u − ū) + (v − v̄)
∂x 0 ∂u 0 ∂v 0
= x̄˙ + Ax̃ + Bũ + Lṽ (2.11)
Se tiene una ecuación lineal para x̃˙ en términos de x̃, ũ y v. Una ecuación lineal
para las desviaciones de los estados y control de sus valores nominales. Siempre
que las desviaciones permanezcan pequeñas, la linealización será exacta y la
ecuación lineal describirá con precisión las desviaciones de x de sus valores no-
minales x̄. En una forma similar se puede expandir la ecuación de medición no
lineal dada por la Ecuación (2.5) alrededor de un punto de operación nominal
x = x̄ y w = w̄ = 0. Esto resulta en la ecuación de medición linealizada
∂h ∂h
ỹ = x̃ + w̃
∂x 0 ∂w 0
= Cx̃ + Dw (2.14)
x̃ = x − x̄; ũ = u − ū; ỹ = y − ȳ
2.4. Discretización 51
2.4. Discretización
Entonces,
donde
∫ T [ ]
−Aτ (AT)2 (AT)3 (AT)4
e dτ = AT − + − + · · · A−1
2! 3! 4!
0
[ ( 2
)]
(AT) (AT)3 (AT)4
= I − I − AT + − + − · · · A−1
2! 3! 4!
Entonces,
∫ T
( )
e−Aτ dτ = I − e−AT A−1 (2.18)
0
F = eAT
( ) (2.19)
G = F I − e−AT A−1 B.
F = eaT
( )b b ( aT )
G = eaT 1 − e−aT = e −1
a a
El sistema de muestreo se convierte ası́ en
b ( aT )
xk = eaT xk−1 + e − 1 uk−1
a
[ ] [1 ]
1 T T2
xk = xk−1 + 2 uk−1
0 1 T
y = [1 0] xk−1
54 Capı́tulo 2. Elementos de sistemas dinámicos
2.5. Estabilidad
En esta sección, se revisará el concepto de estabilidad para sistemas lineales
invariantes en el tiempo. En la primera parte se analizan los sistemas en tiem-
po continuo y luego los sistemas en tiempo discreto. Se establecerán resultados
importantes sin demostración. Para más detalle y resultados adicionales se re-
fiere al lector interesado a algunas referencias estándar sobre sistemas lineales,
v.gr.: [Åström and Murray, 2008], [Chen, 1984], [Jamshidi et al., 1992], [Kailath
et al., 2000].
ẋ = Ax
(2.20)
y = Cx
para alguna matriz constante M. Esto es solo una forma de decir que la matriz
exponencial no aumenta ilimitadamente.
Recordando que eAt = QeÂt Q−1 , donde Q es una matriz constante que contie-
ne los autovectores de A y  es la forma de Jordan de A. El exponencial eÂt
por lo tanto contiene términos como e(λi t) , te(λi t) , t2 e(λi t) , etc, donde λi , es un
autovalor de A. La acotación de eAt es por lo tanto referida a los autovalores
de A como se establece por los siguientes teoremas.
56 Capı́tulo 2. Elementos de sistemas dinámicos
Teorema 2.3
Un sistema lineal de tiempo continuo, invariante en el tiempo es marginalmente
estable si y solo si una de las siguientes condiciones se mantiene.
1. Todos los autovalores de A tienen parte real negativa.
2. Todos los autovalores de A tienen parte real negativa o cero, y esos con
parte real igual a cero tienen una multiplicidad geométrica igual a su multi-
plicidad algebraica. Esto es, los bloques de Jordan que están asociados con
los autovalores que tienen parte real igual a cero son de primer orden.
Teorema 2.4
Un sistema lineal de tiempo continuo, invariante en el tiempo es asintóticamente
estable si y solo si todos los autovalores de A tienen parte real negativa.
Ejemplo 2.3
Considere el sistema
[ ]
0 1 0
ẋ = 0 0 0 x
0 0 −1
Ya que la matriz A es triangular superior, se sabe que sus autovalores están sobre
la diagonal; esto es, los autovalores de A son iguales a 0, 0, y −1. Se ve que el
sistema es asintóticamente inestable puesto que algunos de sus autovalores son no
negativos. También nótese que la matriz A ya está en la forma de Jordan y se
ve que el bloque de Jordan correspondiente al autovalor 0 es de segundo orden.
Por lo tanto, el sistema también es marginalmente inestable. La solución de este
sistema es
x(t) = exp(At)x(0)
[ ]
1 t 0
= 0 1 0 x(0)
0 0 e−t
y x(t) será acotado para todo t. Sin embargo, esto nada dice acerca de la estabi-
lidad del sistema; sólo dice que existe algún x(0) que resulta en un x(t) acotado.
Si se escoge x(0) = [0 1 0]⊤ , entonces
x(t) = exp(At)x(0)
[ ][ ] [ ]
1 t 0 0 t
= 0 1 0 1 = 1
0 0 e−t 0 0
y x(t) crecerá sin lı́mite. Esto prueba que el sistema es asintóticamente inestable
y marginalmente inestable.
Ejemplo 2.4
Considérese el sistema
[ ]
0 0 0
ẋ = 0 0 0 x
0 0 −1
Los autovalores de A son iguales a 0, 0, y −1. Se ve que el sistema es asintótica-
mente inestable puesto que algunos de sus autovalores son no negativos. Además
para ver si el sistema es marginalmente estable, se necesita calcular la multipli-
cidad geométrica del autovalor 0. (Esto puede hacerse notando que A ya está en
la forma de Jordan, pero se irá a través del ejercicio por motivo de ilustración).
resolviendo la ecuación
[ ]
0
(λI − A)v = 0
0
y x(t) se aproxima a 0 cuando t crece. Sin embargo, esto no dice nada acerca de la
estabilidad asintótica del sistema; esto solamente dice que existe algún x(0) que
resulta en el estado x(t) que se aproxima asintóticamente a 0. Si por el contrario
se opta por escoger x(0) = [0 1 0]⊤ , entonces
[ ][ ] [ ]
1 0 0 0 0
x(t) = 0 1 0 1 = 1
0 0 e−t 0 0
xk+1 = Fxk
(2.25)
yy = Hxk
estable si el estado xk es acotado para todo k y para todos los estados iniciales
acotados x0 .
Definición 2.4
Un sistema lineal de tiempo discreto invariante en el tiempo es asintóticamente
estable si,
lı́m xk = 0 (2.26)
k→∞
xk = Fk x0 (2.27)
para alguna matriz constante M. Esto es sólo una forma de decir que las potencias
de F no aumenta sin lı́mite.
Teorema 2.6
Un sistema lineal de tiempo discreto, invariante en el tiempo es asintóticamente
estable si y solo si
lı́m Fk = 0 (2.29)
k→∞
Recordando que Fk = QF̂k Q−1 , donde Q es una matriz constante que contiene
los autovectores de F y F̂ es la forma de Jordan de F. La matriz F̂k por lo tanto
contiene términos como λki , kλki , k 2 λki , etc, donde λi , es un autovalor de F. La
60 Capı́tulo 2. Elementos de sistemas dinámicos
Teorema 2.8
Un sistema lineal de tiempo discreto, invariante en el tiempo es asintóticamente
estable si y solo si todos los autovalores de F tienen magnitud menor que uno.
2.6.1. Controlabilidad
Las siguientes definiciones y teoremas dan descripción clara de la controlabili-
dad para sistemas lineales tanto en tiempo continuo como en discreto.
Definición 2.5
Un sistema de tiempo continuo es controlable si para cualquier estado inicial x(0)
2.6. Controlabilidad y Observabilidad 61
y cualquier tiempo final t > 0 existe una ley de control que transfiera el estado a
cualquier valor deseado en el tiempo t.
Definición 2.6
Un sistema de tiempo discreto es controlable si para cualquier estado inicial x0 y
algún tiempo final k existe una ley de control que transfiera el estado a cualquier
valor deseado en el tiempo k.
Teorema 2.10
Un sistema lineal de n estados, continuo e invariante en el tiempo ẋ = Ax + Bu
es controlable si y solo si la matriz grammiana de controlabilidad definida por
∫ t
⊤
G = eAτ BB⊤ eA τ dτ (2.31)
0
Teorema 2.11
Un sistema lineal de n estados, continuo e invariante en el tiempo ẋ = Ax + Bu
es controlable si y solo si la ecuación diferencial de Lyapunov
tiene una solución definida positiva W(t) para algún t ∈ (0, ∞). Ésta también es
llamada ecuación de Sylvester.
Teorema 2.13
Un sistema lineal de n estados, discreto e invariante en el tiempo xk = Fxk−1 +
Guk−1 es controlable si y solo si la matriz grammiana controlabilidad definida
por
∑
k
G = Fk−1 GG⊤ (A⊤ )k−i (2.34)
i=0
Teorema 2.14
Un sistema lineal de n estados, discreto e invariante en el tiempo xk = Fxk−1 +
Guk−1 es controlable si y solo si la ecuación de diferencias de Lyapunov
tiene una solución definida positiva Wk para algún k ∈ (0, ∞). Ésta es también
llamada la ecuación de Stein.
Nótese que los Teoremas 2.9 y 2.12 dan pruebas idénticas de controlabilidad
tanto para sistemas de tiempo continuo como de tiempo discreto. En general,
estas son las pruebas más simples de controlabilidad. Las pruebas de contro-
labilidad para sistemas lineales variantes en el tiempo pueden obtenerse gene-
ralizando los teoremas anteriores. La controlabilidad para sistemas no lineales
es mucho más difı́cil de formalizar.
2.6. Controlabilidad y Observabilidad 63
Ejemplo 2.5
El circuito RLC de la Figura 2.1 tiene la descripción en variables de estados dada
por (ver por ejemplo [Chan et al., 1972])
[ ] [ −2 1 ] [ ] [ 1 ]
v̇C RC C
vC RC
= −1 + 1 u
ı̇L L 0 ıL L
+ R
u L R
-
2.6.2. Observabilidad
Las siguientes definiciones y teoremas dan descripción clara de la observabili-
dad para sistemas lineales tanto en tiempo continuo como en discreto.
Definición 2.7
Un sistema de tiempo continuo es observable si para cualquier estado inicial x(0) y
cualquier tiempo final t > 0 el estado inicial x(0) se puede determinar únicamente
con el conocimiento de la entrada u(τ ) y la salida y(τ ) para todo τ ∈ [0, t].
Definición 2.8
Un sistema de tiempo discreto es observable si para cualquier estado inicial x0 y
algún tiempo final k el estado inicial x0 se puede determinar únicamente con el
conocimiento de la entrada ui y la salida yi para todo i ∈ [0, k].
Teorema 2.16
El sistema lineal de n estados, continuo e invariante en el tiempo
ẋ = Ax + Bu
(2.39)
y = Cx
es observable si y solo si la observabilidad grammiana definida por
∫ t
⊤
eA τ C⊤ CeAτ dτ (2.40)
0
−Ẇ = WA + A⊤ W + C⊤ C (2.42)
W(t) = 0
tiene una solución definida positiva W(τ ) para algún τ ∈ (0, t). Esta es también
llamada la ecuación de Sylvester.
Similar al caso de tiempo continuo, los siguientes teoremas pueden usarse para
probar la observabilidad de un sistema lineal discreto invariante en el tiempo.
Teorema 2.18
El sistema lineal de n estados, discreto e invariante en el tiempo
xk = Fxk−1 + Guk−1
(2.43)
yk = Hxk
tiene la matriz de observabilidad O defina por
H
HF
O= ...
(2.44)
HFn−1
Teorema 2.19
Un sistema lineal de n estados, discreto e invariante en el tiempo
xk = Fxk−1 + Guk−1
(2.45)
yk = Hxk
es observable si y solo si la observabilidad grammiana definida por
∑
k
(F⊤ )i H⊤ HFi (2.46)
i=0
Teorema 2.20
El sistema lineal de n estados, discreto e invariante en el tiempo
xk = Fxk−1 + Guk−1
(2.47)
yk = Hxk
es observable si y solo si la ecuación diferencial de Lyapunov
tiene una solución definida positiva W0 para algún k ∈ (0, ∞). Esta es también
llamada la ecuación de Stein.
Nótese que los Teoremas 2.15 y 2.18 dan pruebas idénticas de observabilidad
tanto para sistemas de tiempo continuo como de tiempo discreto. En general,
estas son las pruebas más simples de observabilidad. Las pruebas de observa-
bilidad para sistemas lineales variantes en el tiempo pueden obtenerse genera-
lizando los anteriores teoremas. La observabilidad para sistemas no lineales es
mucho más difı́cil de formalizar.
Ejemplo 2.6
El circuito RLC de la Figura 2.1 está definido en variables de estados como
[ ] [ −2 1 ] [ ] [ 1 ]
v̇C RC C
vC RC
= −1 + 1 u
ı̇L 0 ıL
L
[ ] L
vC
y = [−1 0]
ıL
2.6. Controlabilidad y Observabilidad 67
Por lo tanto el sistema es no observable. Serı́a muy difı́cil obtener este resultado
a partir de los Teoremas 2.16 y 2.17.
ẋ = Ax + Bu
(2.49)
y = Cx + Du
ẋ = Ax + Bu
[ ] [ ]
1 1 2 1
= 0 1 3 x+ 1 u
0 0 −2 0
y = Cx + Du
= [1 0 0] x
Este sistema tiene la misma función de transferencia que
M = [v1 v2 v3 ]
[ ]
1 0 1
= 0 1 3
0 0 −3
Se pueden observar las equivalencias:
Ā = M−1 AM
B̄ = M−1 B
C̄ = CM
El sistema de la forma de Jordan tiene dos modos desacoplados.
El primer modo es
[ ] [ ]
1 1 1
x̄˙ = x̄ + u
0 1 1 1
y1 = [1 0] x̄1
El segundo modo es
x̄˙ 2 = −2x̄2 + 0u
(2.52)
y2 = x̄2
Definición 2.9
Si un sistema es controlable o estable, entonces también es estabilizable. Si un
sistema es no controlable o inestable, entonces es estabilizable si sus modos no
controlables son estables.
Definición 2.10
Si un sistema es observable o estable, entonces también es detectable. Si un sistema
es no observable o inestable, entonces es detectable si sus modos no observables
son estables.
70 Capı́tulo 2. Elementos de sistemas dinámicos
2.7. Resumen
En este capı́tulo se ha hecho una revisión de aspectos básicos de álgebra lineal
incluyendo una introducción a productos de Hadamard y de Kronecker, los
cuales no son muy corrientes en los cursos introductorios. La razón es que se
presentan algunas procesos en la práctica que requieren esas aplicaciones, tales
como en la solución de la ecuación de Sylvester o Lyapunov que tiene ocurrencia
en sistemas de control. Esto conduce a un conjunto de resultados que se han
expuesto en el texto en forma de teoremas con el objeto de darles formalidad
y también como ejercicio para el lector. También se ha hecho un estudio del
concepto de número de condición. Esto es de utilidad en el proceso de análisis
de estabilidad numérica de algunos sistemas dinámicos. En la segunda parte
se abordó el problema del cálculo de matrices, procedimiento que conduce a
la obtención de resultados directos para diversas aplicaciones de matrices y
vectores. Se ha procurado realizar los procedimientos por el método de paso
a paso con el fin de que el lector pueda seguir fácilmente las demostraciones.
Finalmente, se hace una pequeña introducción a los sistemas lineales y no
lineales y a la discretización de los procesos dinámicos.
Bibliografı́a
Chan, S.-P., Chan, S.-Y., and Chan, S.-G. (1972). Analysis of Linear Networks and Systems.
A Matrix-Oriented Approach with Computer Applications. Addison-Wesley Series in
ELectrical Engineering. Addison-Wesley Publishing Company, Reading, MA.
Chen, C.-T. (1984). Linear System Theory and Design. Saunders College Pu., Harcourt
Brace Jovanovich, Inc., Orlando, FL 32887, USA, 2a edition.
Jamshidi, M., Tarokh, M., and Shafai, B. (1992). Computer-Aided Analysis and Design
of Linear Control Systems. Prentice-Hall International, Inc., A Simon and Schuster
Company Englewood Cliffs, N J 07632, USA, 1a edition.
Kailath, T., Sayed, A., and Hassibi, B. (2000). Linear Estimation. Prentice-Hall, Upper
Saddle River, New Jersey, USA.
Kalman, R. E. (1959). On the general theory of control systems. IRE Transactions on
Automatic Control, 4(3): 110–122.
Kalman, R. E. (1960). Contributions to the theory of optimal control. Boletı́n de la Sociedad
Matemática Mexicana, 5: 102–119.
Khalil, H. (1996). Nonlinear Systems. Simon and Schuster, Viacom Company Upper Saddle
River, N. J. 07458, USA, segunda edition.
Lewis, F. L. (1992). Applied Optimal Control & Estimation: Digital Design & Implemen-
tation. Digital Signal Processing. Prentices-Hall, Inc., Simon & Schuster Company
Englewoods Cliffs, NJ 07632, USA.
Simon, D. (2006). Optimal State Estimation: Kalman, H-infinity, and Nonlinear Approaches.
John Wiley & Sons, Inc., New Jersey, USA.
Åström, K. J. and Murray, R. (2008). Feedback systems: an introduction for scientists and
engineers. Princeton University Press, 41 William Street, Princeton, New Jersey 08540,
USA.
Åström, K. J. and Wittenmark, B. (1997). Computer-Controlled Systems: Theory and De-
sign. Prentice Hall Information and System Sciences, 3a edition.
71
Capı́tulo 3
3.1. Introducción
73
74 Capı́tulo 3. Estimación por mı́nimos cuadrados
y1 = h11 x1 + · · · + h1n xn + w1
.. ..
. .
yk = hk1 x1 + · · · + hkn xn + wk
y = Hx + w (3.1)
ϵy = y − Hx̂ (3.2)
J = ϵ2y1 + · · · + ϵ2yk
= ϵ⊤
y ϵy (3.3)
3.2. Estimación de una constante 75
y⊤ H = x̂⊤ H⊤ H
H⊤ y = H⊤ Hx̂
Es decir,
∂ 2J
= 2H⊤ H (3.6)
∂x̂2
Si H es una matriz de orden k × n, entonces se verifica que tanto HH⊤ como
H⊤ H son matrices semidefinidas positivas. Ahora bien,
Si ρ(H) = k entonces HH⊤ es definida positiva.
Si ρ(H) = n entonces H⊤ H es definida positiva.
Este último caso corresponde a la Ecuación (3.6) lo cual completa la demos-
tración.
76 Capı́tulo 3. Estimación por mı́nimos cuadrados
Ejemplo 3.1
La idea es tratar de estimar la resistencia x de un resistor no marcado con base a
k medidas ruidosas de un multı́metro. En este caso x es un escalar de modo que
las k medidas ruidosas estarán dadas como
y 1 = x + w1
.. ..
. .
y k = x + wk
Estas k ecuaciones se pueden combinar en una simple expresión vectorial como
y1 1
.. = ... x + w
.
yk 1
1∑
k
x̂ = yi (3.7)
k
i=1
En este simple ejemplo, se puede observar que la estimación por mı́nimos cua-
drados concuerda con nuestra intuición de simplemente calcular el promedio
de las medidas.
E {wi2 } = σi2 i = 1, . . . , k
ϵ2y1 ϵ2yk
J= + · · · + (3.8)
σ12 σk2
J = ϵ⊤ −1 ⊤ −1
y R ϵy = (y − Hx̂) R (y − Hx̂)
es decir,
Finalmente,
y1 1 w1
... = ... x + ...
yk 1 wk
R = diag(σ12 · · · σk2 )
3.4. Estimación recursiva por mı́nimos cuadrados 79
yk = H k x + w k
(3.10)
x̂k = x̂k−1 + Kk (yk − Hk x̂k−1 )
Para esto, se calcula x̂ sobre la base del valor estimado anterior x̂k−1 y la
nueva medida yk . Kk es una matriz a ser determinada, llamada la matriz
de ganancia del estimador. La cantidad (yk − Hk x̂k−1 ), se denomina término
corrector. Nótese que si el término corrector es cero o si la matriz de ganancia
es cero, entonces el valor estimado no cambia del paso (k − 1) al paso k. Antes
de calcular el valor óptimo de la matriz de ganancia Kk , se debe analizar el
promedio del error de estimación del estimador recursivo lineal.
El promedio del error estimado se puede calcular como
E {ϵx,k } = E {x − x̂k }
= E {x − x̂k−1 − Kk (yk − Hk x̂k−1 )}
= E {ϵx,k−1 − Kk (Hk x + wk − Hk x̂k−1 )}
= E {ϵx,k−1 − Kk Hk (x − x̂k−1 ) − Kk wk }
= (I − Kk Hk )E {ϵx,k−1 } − Kk E {wk } (3.11)
80 Capı́tulo 3. Estimación por mı́nimos cuadrados
Jk = E {ϵ⊤ ⊤ ⊤
x,k ϵx,k } = E {Tr(ϵx,k ϵx,k )} = Tr(E {ϵx,k ϵx,k })
= Tr(Pk ) (3.13)
Pk = E {ϵx,k ϵ⊤
x,k }
puesto que ambos valores esperados son cero. Por lo tanto, la Ecuación (3.14)
se reduce a
Pk = (I − Kk Hk )Pk−1 (I − Kk Hk )⊤ + Kk Rk K⊤
k (3.16)
∂Jk ∂Tr(Pk )
= = 2(I − Kk Hk )Pk−1 (−H⊤
k ) + 2Kk Rk (3.17)
∂Kk ∂Kk
Ahora bien, para encontrar el valor de Kk que minimice Jk , se hace la derivada
igual a cero y entonces se resuelve para Kk como sigue:
Kk Rk = (I − Kk Hk )Pk−1 H⊤
k
Kk Rk = Pk−1 H⊤ ⊤
k − Kk Hk Pk−1 Hk
82 Capı́tulo 3. Estimación por mı́nimos cuadrados
Kk (Rk + Hk Pk−1 H⊤ ⊤
k ) = Pk−1 Hk
Kk = Pk−1 H⊤ ⊤
k (Hk Pk−1 Hk + Rk )
−1
(3.18)
Las ecuaciones (3.10), (3.16), y (3.18) forman la estimación recursiva por mı́-
nimos cuadrados, la cual se resume en el Algoritmo 3.1.
for k = 1 to N do
yk = Hk xk + wk ◃ wk ∼ N ID(0, Rk )
end for
2. Actualización del valor estimado de x y la covarianza del error de estimación P.
for k = 1 to N do
( )−1
Kk = Pk−1 H⊤ ⊤
k Hk Pk−1 Hk + Rk
x̂k = x̂k−1 + Kk (yk − Hk x̂k−1 )
Pk = (I − Kk Hk )Pk−1 (I − Kk Hk )⊤ + Kk Rk K⊤
k
end for
Pk = Pk−1 − 2Pk−1 H⊤ −1 ⊤ −1 ⊤ −1
k Sk Hk Pk−1 + Pk−1 Hk Sk · Hk Pk−1 Hk · Sk Hk Pk−1
+ Pk−1 H⊤ −1 −1
k Sk · Rk · Sk Hk Pk−1 .
−1 ⊤ −1 −1
Pk = Pk−1 − 2Pk−1 H⊤ ⊤
k Sk Hk Pk−1 + Pk−1 Hk Sk (Hk Pk−1 Hk + Rk )Sk Hk Pk−1
= Pk−1 − 2Pk−1 H⊤ −1 ⊤ −1
k Sk Hk Pk−1 + Pk−1 Hk Sk Hk Pk−1
−1
= Pk−1 − Pk−1 H⊤
k Sk Hk Pk−1 (3.20)
Pk = Pk−1 − Kk Hk Pk−1
= (I − Kk Hk )Pk−1 (3.21)
Comparada con la Ecuación (3.16), esta es una ecuación más simple para la
determinación de Pk ; sin embargo, esta expresión puede causar problemas de
cálculo numérico (v.gr.: de tipo escalamiento), por no ser definida positiva, aun
cuando Pk−1 y Rk sı́ lo sean.
Otra forma alterna consiste en la aplicación del lema de inversión de matrices
(Ecuación (1.72)), para reescribir la ecuación de la medida actualizada de Pk .
Comenzando con la Ecuación (3.20) se obtiene:
Pk = Pk−1 − Pk−1 H⊤ ⊤ −1
k (Hk Pk−1 Hk + Rk ) Hk Pk−1 (3.22)
Pk = Pk−1 − Pk−1 H⊤ ⊤ −1
k (Hk Pk−1 Hk + Rk ) Hk Pk−1
= A−1 − A−1 B(CA−1 B + D)−1 CA−1
= (A + BD−1 C)−1 ≡ (P−1 ⊤ −1 −1
k−1 + Hk Rk Hk ) , (3.23)
es decir,
( )−1
Pk = P−1 ⊤ −1
k−1 + Hk Rk Hk (3.24)
Esta ecuación es más complicada para el cálculo de Pk puesto que requiere tres
inversiones de matrices, pero puede ser más ventajosa en algunos casos, como
se verá más adelante. Se puede utilizar la Ecuación (3.24) con el fin de obtener
una ecuación equivalente para la ganancia del estimador Kk . Reescribiendo la
Ecuación (3.18):
Kk = Pk−1 H⊤ ⊤
k (Hk Pk−1 Hk + Rk )
−1
(3.25)
Kk = Pk P−1 ⊤ ⊤
k Pk−1 Hk (Hk Pk−1 Hk + Rk )
−1
(3.26)
Sustituyendo P−1
k a partir de la Ecuación (3.24):
Kk = Pk (P−1 ⊤ −1 ⊤ ⊤
k−1 + Hk Rk Hk )Pk−1 Hk (Hk Pk−1 Hk + Rk )
−1
(3.27)
Kk = Pk (H⊤ ⊤ −1 ⊤ ⊤
k + Hk Rk Hk Pk−1 Hk )(Hk Pk−1 Hk + Rk )
−1
(3.31)
Ahora se extrae H⊤
k del paréntesis, con lo cual se tiene:
Kk = Pk H⊤ −1 ⊤ ⊤
k (I + Rk Hk Pk−1 Hk )(Hk Pk−1 Hk + Rk )
−1
(3.32)
3.5. Formas alternas del estimador 85
x̂0 = E {x}
(3.29)
P0 = E {(x − x̂0 )(x − x̂0 )⊤ }
Para k ∈ {1, 2, · · · , ∞}, las ecuaciones de estimación por mı́nimos cuadrados
son:
Kk = Pk−1 H⊤ ⊤
k (Hk Pk−1 Hk + Rk )
−1
= Pk H ⊤ −1
k Rk
x̂k = x̂k−1 + Kk (yk − Hk x̂k−1 )
(3.30)
Pk = (I − Kk Hk )Pk−1 (I − Kk Hk )⊤ + Kk Rk K⊤
k
= (P−1 ⊤ −1
k−1 + Hk Rk Hk )
−1
= (I − Kk Hk )Pk−1
Kk = Pk H⊤ −1 ⊤ ⊤
k Rk (Rk + Hk Pk−1 Hk )(Hk Pk−1 Hk + Rk )
−1
Kk = Pk H⊤ −1
k Rk (3.33)
yk = Hk x + wk
Hk = 1 (3.34)
Rk = E {wk2 }
Para este problema, la matriz de medida Hk es un escalar y la covarianza del
ruido de la medida Rk también es un escalar. Se supone que cada medida tiene la
misma covarianza, ası́ la covarianza del ruido de la medida Rk no es una función
de k y puede escribirse como R. Inicialmente, antes de tomar cualquier medida,
se tiene alguna idea del valor de la resistencia x y esto constituye el estimado
inicial. También se tiene alguna incertidumbre acerca del estimado inicial y esto
conforma la covarianza inicial:
x̂0 = E {x}
P0 = E {(x − x̂0 )(x − x̂0 )⊤ } (3.35)
P0 = E {(x − x̂0 )2 }
Si no se tiene idea en absoluto del valor de la resistencia, entonces P0 → ∞. Si
se está 100 % seguro acerca del valor de la resistencia antes de tomar cualquier
medida, entonces P0 = 0. El conjunto de ecuaciones (3.30) dice como obtener la
ganancia del estimador, el estimado de x y la covarianza de la estimación, después
de la primera medición (k = 1):
P1 P0
K2 = =
P1 + R 2P0 + R
P1 P0 + R P0
x̂2 = x̂1 + (y2 − x̂1 ) = x̂1 + y2 (3.36)
P1 + R 2P0 + R 2P0 + R
P1 R P0 R
P2 = =
P1 + R 2P0 + R
3.5. Formas alternas del estimador 87
P0
Kk = (3.37)
kP0 + R
(k − 1)P0 + R P0 k−1 1
x̂k = lı́m x̂k−1 + yk = x̂k−1 + yk
P0 →∞ kP0 + R kP0 + R k k
1
= [(k − 1)x̂k−1 + yk ] (3.39)
k
En otras palabras, el valor estimado óptimo de x es igual al promedio de las
medidas yk , lo cual se puede escribir como
1 ∑ k
1 ∑
k−1
1 1 ∑
k−1
ȳk = yj = yj + yk = (k − 1) yj + yk
k k k k−1
j=1 j=1 j=1
1
= [(k − 1)ȳk−1 + yk ] (3.40)
k
Ejemplo 3.4
Supóngase que se desea ajustar una lı́nea recta a un conjunto de datos puntua-
les. El problema de ajuste datos lineales se puede escribir como
yk = x1 + x2 tk + wk
(3.41)
E {wk2 } = Rk
Kk = Pk−1 H⊤ ⊤
k (Hk Pk−1 Hk + Rk )
−1
Ejemplo 3.5
Supóngase que un tanque contiene una concentración x1 del quı́mico 1 y una
concentración x2 del quı́mico 2. Se tiene una instrumentación que puede detectar
la concentración combinada (x1 +x2 ) de los dos quı́micos, pero la instrumentación
no puede distinguir ninguno de los quı́micos por separado.
Se elimina del tanque el quı́mico 2 a través de un proceso de lixiviación de
modo que su concentración decrece en un 1 % por cada medida. La ecuación de
medición está dada por
yk = x1 + 0.99k−1 x2 + wk
[ ]
= 1 0.99k−1 x + wk (3.44)
11
10
9
Estimados
8 x1
7 x2
6
5
4
0 10 20 30 40 50 60
1
P(1,1)
0.8 P(2,2)
Varianza
0.6
0.4
0.2
0
0 10 20 30 40 50 60
Paso de tiempo
como
∑
∞
SXY (ω) = RX (k)e−jωk ω ∈ [−π, π]
−∞
∫ ∞
1
RXY (ω) = SX (ω)ejkω dω (3.48)
2π −∞
RX (ω) =σ 2 δk (3.49)
RX (τ ) = RX (0)δ(τ ) (3.52)
Para encontrar el filtro óptimo G(ω) se necesita minimizar E {e2 (t)}, lo cual
significa que se necesita conocer Sx (ω) y Sv (ω), las propiedades estadı́sticas de
la señal x(t) y el ruido v(t).
de las entradas presentes y futuras. Los sistemas del mundo real son siempre
causales, pero un filtro que se utiliza para postprocesado puede ser no causal.)
1
G(ω) = (3.67)
1 + T jω
Esto puede ser una suposición no válida, pero reduce el problema a un problema
de optimización paramétrica. Con el fin de simplificar aún más el problema,
supóngase que Sx (ω) y Sv (ω) están dados de la siguiente manera:
2σ 2 β
Sx (ω) =
ω2 + β 2 (3.68)
Sv (ω) = A
Ahora se toma un enfoque más general para encontrar el filtro óptimo. El valor
esperado del error de estimación puede calcularse como
∫ ∫
+ g(u)g(γ)[x(t − u) + v(t − u)][x(t − v) + v(t − v)]dudγ
(3.71)
∫ ∫
+ g(u)g(γ)[Rx (u − v) + Rv (u − v)]dudγ (3.72)
Ahora se puede usar un enfoque del cálculo variacional [Weinstock, 1974] para
encontrar el filtro g(t) que minimize E {e2 (t)}. Se reemplaza g(t) en la anterior
ecuación con g(t) + ϵη(t), donde ϵ es algún número pequeño, y η(t) es una per-
turbación arbitraria en g(t). Aplicando los conceptos del cálculo de variaciones
los cuales establecen que se puede minimizar E {e2 (t)} haciendo
∂E {e2 (t)}
=0 (3.73)
∂ϵ ϵ=0
3.6. El filtro de Wiener 97
Esto da la condición necesaria para la optimalidad del filtro g(t) como sigue:
∫ [ ∫ ]
η(τ ) −Rx (τ ) + g(u)[Rx (u − τ ) + Rv (u − τ )]du dτ = 0 (3.77)
Se requiere resolver esta ecuación para g(t) y ası́ encontrar el filtro óptimo.
98 Capı́tulo 3. Estimación por mı́nimos cuadrados
∫
Rx (τ ) = g(u)[Rx (u − τ ) + Rv (u − τ )]du = g(τ ) ∗ [Rx (τ ) + Rv (τ )]
2
Sx (ω) =
ω2 +1
Sv (ω) = 1
( √ )
2 1 2 3
G(ω) = 2 =√
ω +3 3 ω2 + 3
1 √
g(t) = √ e− 3|t| ≈ 0.58e−0.58|t| , t ∈ (−∞, ∞)
3
Además, para obtener una representación en el dominio del tiempo del filtro, se
realiza una expansión en fracciones parciales de G(ω) a fin de encontrar la parte
causal y la parte anti-causal del filtro:
3.6. El filtro de Wiener 99
1 1
G(ω) = √ √ +√ √
3(jω + 3) 3(−jω + 3)
| {z } | {z }
filtro causal filtro anti-causal
De esto se ve que
1 1
X̂(ω) = √ √ Y (s) − √ √ Y (s)
3(jω + 3) 3(jω − 3)
= X̂c (ω) + X̂a (ω)
X̂c (ω) y X̂a (ω) (definidos por la anterior ecuación) son la parte causal y anti-
causal de X̂(ω), respectivamente. En el dominio del tiempo, esto puede escribirse
como
x̂(t) = x̂c (t) + x̂a (t)
√ √
x̂˙ c (t) = − 3x̂c + y/ 3
√ √
x̂˙ a (t) = 3x̂a − y/ 3
La ecuación x̂˙ c (t) corre hacia adelante en el tiempo y es por lo tanto causal y
estable. La ecuación x̂˙ a (t) corre hacia atrás en el tiempo y es por lo tanto anti-
causal y estable. (Si ésta corriera hacia adelante en el tiempo, serı́a inestable).
Sx (ω) A(ω)
+
G(ω)Sxv (ω) = −
− − (3.83)
Sxv (ω) Sxv (ω)
+ Sx (ω)
G(ω)Sxv (ω) = parte causal de −
Sxv (ω)
[ ] (3.84)
1 Sx (ω)
G(ω) = + (ω)
parte causal de − (ω)
Sxv Sxv
Ejemplo 3.9
Considérese el sistema discutido en el Ejemplo 3.8, con A = β = σ = 1.
2
Sx (ω) =
ω2 +1
Sv (ω) = 1
ω2 + 3
Sxv (ω) =
ω2 + 1
Dividiendo esto en sus factores causal y anti-causal da
( √ )( √ )
jω + 3 −jω + 3
Sxv (ω) =
jω + 1 −jω + 1
| {z }| {z }
+ −
Sxv (ω) Sxv (ω)
Sx (ω) 2(−jω + 1)
− = √
Sxv (ω) 2
(ω + 1)(−jω + 3)
2
= √
(−jω + 3)(jω + 1)
√ √
3−1 3−1
= + √
jω + 1 −jω + 3
| {z } | {z }
parte causal parte anti-causal
La Ecuación (3.84) da
( ) (√ ) ( √ )
jω + 1 3−1 3−1
G(ω) = √ = √
jω + 3 jω + 1 jω + 3
√ √
g(t) = ( 3 − 1)e− 3t , t≥0
3.6.5. Comparación
Comparando los tres ejemplos de diseño de filtro óptimo presentados en esta
sección (Ejemplos 3.7, 3.8 y 3.9), se puede demostrar que los valores esperados
de los errores cuadráticos medios de la filtro son las siguientes [Brown and
Hwang, 2012]:
Método de optimización de parámetros: E {e2 (t)} = 0.914
102 Capı́tulo 3. Estimación por mı́nimos cuadrados
3.7. Resumen
En este capı́tulo se discutió la estimación por mı́nimos cuadrados en un par de
contextos diferentes. Primero se derivó un método para estimar vector cons-
tante sobre la base de varias mediciones ruidosas de ese vector. De hecho, las
mediciones no tienen que ser medidas directas del vector constante, pero pue-
den ser mediciones de alguna combinación lineal de los elementos de ese vector
constante. Además, el ruido asociado con cada medición no tiene que ser el
mismo. En la técnica de estimación por mı́nimos cuadrados que se deriva, se
asume que el ruido de medición es blanco de media cero (no correlacionado de
un paso de tiempo al siguiente) y que se conoce su varianza. A continuación,
se extiende el estimador de mı́nimos cuadrados a una formulación recursiva,
en la que el esfuerzo de cálculo sigue siendo el mismo en cada paso de tiempo,
independientemente del número total de mediciones que sean procesadas. La
estimación por mı́nimos cuadrados de un vector constante forma parte de la
fundamentación del desarrollo del filtro de Kalman, que se va a derivar más
adelante.
Finalmente, se da un breve resumen de los conceptos correspondientes a los
filtros de Wiener, que si bien plantearon paradigmas difı́ciles de superar, hoy
dı́a no se aplican ampliamente debido a dificultades en el desarrollo de los
procedimientos, ya que para su aplicación se requieren señales estacionarias
y eventualmente monovariables. Para los casos donde hay señales estadı́sticas
variables en el tiempo o de variables múltiples, el mejor estimador es el filtro
de Kalman.
Bibliografı́a
103
Capı́tulo 4
4.1. Introducción
105
106 Capı́tulo 4. Filtro de Kalman discreto
1. Existe una solución única para P si y solo si λi (F)λj (F) ̸= 1 para todo i, j.
Esta solución única es simétrica.
∑
k−1
xk = Fk,0 x0 + (Fk,i+1 vi + Fk,i+1 Gi ui ) (4.6)
i=0
xk ∼ N (x̄k , Pk ). (4.8)
Ejemplo 4.1
Considérese el sistema
[1 ]
2 0
xk+1 = x + vk
0 12
vk ∼ (0, Q)
[ ]
1 1
Q=
0 1
Solución:
∞
∑ ∞ [1
∑ ][ ]i [ ]i
i ⊤ i 1 1 12 0
2 0
P= F Q(F ) = =
0 12 0 1 0 1
2
i=0 i=0
[ ] [1 1] [ 1 1
] [1 1
]
1 1
1 + ···
= + 4 4 + 16 16 64 64
1 +
0 1 0 14 0 16 0 64
[ ] ( )0 [ ] ( )1 [ ] ( )2 [ ] ( )3
1 1 1 1 1 1 1 1 1 1 1 1
= + + + + ···
0 1 4 0 1 4 0 1 4 0 1 4
[ ] ∞ ( )i
1 1 ∑ 1
=
0 1 4
i=0
( )i
1 − ri 1 − 14 4
Si = lı́m = lı́m =
i→∞ 1 − r i→∞ 1 − 1 3
4
4 4
3 3
P=
4
0 3
4.2. Sistemas discretos en el tiempo 109
Ejemplo 4.2
Un sistema lineal describe la población de un predador x1 y la de su presa x2 los
cuales pueden ser expresados como
[ ] [ ]
0.2 0.4 0
xk+1 = xk + u + vk
−0.4 1 1 k (4.10)
vk ∼ (0, Q) Q = diag(1, 2)
20
Población media
x
1
x
2
10
0
0 5 10 15 20 25
Varianza de la población
15 P(1,1)
P(1,2)
10 P(2,2)
0
0 5 10 15 20 25
Pasos de tiempo
Nótese que, puesto que F para este ejemplo es estable y Q es definida positi-
va, entonces el Teorema 4.1 garantiza que P tiene una solución única de estado
estacionario definida positiva.
xk = Fk−1 xk−1 + Gk−1 uk−1 + Lk−1 v̌k−1 , v̌k ∼ (0, Q̌k ) (4.11)
yk = Hk xk + wk , wk ∼ (0, Mk Řk M⊤
k) (4.15)
ẋ = Ax + Bu + v (4.16)
Se sabe que la solución para x(t) en algún tiempo arbitrario tk , está dada por
∫ tk
A(tk −tk−1 )
x(tk ) = e x(tk−1 ) + eA(tk −τ ) [B(τ )u(τ ) + v(τ )]dτ (4.17)
tk−1
Se supone que u(t) = uk para t ∈ [tk , tk+1 ]; esto es, el control u(t) es constante
a tramos.
Definiendo
T = tk − tk−1
xk = x(tk ) (4.18)
uk = u(tk )
∫ tk ∫ tk
A(tk −τ )
xk = e AT
xk−1 + e B(τ )dτ uk−1 + eA(tk −τ ) v(τ )dτ (4.19)
tk−1 tk−1
Fk = eAT
∫ tk+1
(4.20)
Gk = eA(tk+1 −τ ) B(τ )dτ
tk
112 Capı́tulo 4. Filtro de Kalman discreto
Ahora, si se asume que v(t) es ruido blanco de tiempo continuo con una cova-
rianza Qc (t), se ve que
{ }
E v(τ )v⊤ (α) = Qc (τ )δ(τ − α) (4.24)
= Fk−1 Pk−1 F⊤
k−1 + Qk−1 (4.25)
4.3. Sistemas con datos muestreados 113
En general, es difı́cil calcular Qk−1 , pero para valores pequeños de (tk − tk−1 )
se puede obtener
Ejemplo 4.3
Supóngase que se tiene un sistema dinámico de tiempo continuo de primer orden
dado por la ecuación
ẋ = ax + v
(4.28)
E {v(t)v(t + τ )} = qc δ(τ )
Para describir la mayorı́a de los procesos fı́sicos simples se puede emplear ecuacio-
nes de primer orden. Por ejemplo, esta ecuación describe el comportamiento de la
corriente a través de un circuito serie RL que está excitado por una tensión aleato-
ria v(t), donde a = −R/L. Supóngase que se está interesado en obtener la media
y la covarianza del estado x(t) cada T unidades de tiempo; esto es, tk − tk−1 = T.
Para este simple ejemplo escalar, se puede calcular explı́citamente Qk−1 en la
Ecuación (4.26) como
∫ tk
Qk−1 = ea(tk −τ ) qc ea(tk −τ ) dτ
tk−1
∫ tk
=e 2atk
qc e−2aτ dτ
tk−1
[ −2atk−1 ]
2atk e − e−2atk
=e qc
2a
qc [ 2a(tk −tk−1 ) ]
= e −1
2a
qc [ 2aT ]
Qk−1 = e −1 (4.29)
2a
Para valores pequeños de T, se puede expandir la ecuación anterior en una serie
114 Capı́tulo 4. Filtro de Kalman discreto
Ejemplo 4.4
Dos mezclas quı́micas se vierten en un tanque. Una tiene concentración c1 y se
vierte a una razón de F1 y la otra tiene concentración c2 y se vierte a una razón
de F2 . El tanque tiene un volumen V y su salida está a una concentración c y
razón F . Esto es tı́pico de muchos sistemas de control de procesos [Kwakernaak
4.3. Sistemas con datos muestreados 115
and Sivan, 1972]. La ecuación linealizada para este sistema se puede escribir como
[ ] [ ]
− 2V
F0
0 1 1
ẋ = 0 x + c1 −c0 c2 −c0 u + v (4.30)
0 − FV00 V0 V0
Solución:
∫ tk
⊤ (t −τ )
Qk−1 = eA(tk −τ ) Qc eA k
dτ
tk−1
∫ tk
⊤τ ⊤t
=e Atk
· e−Aτ Qc e−A dτ · eA k
(4.31)
tk−1
∫ tk
⊤ (t −τ )
Qk−1 = eA(tk −τ ) Qc eA k
dτ
tk−1
[ ] ∫ tk [ τ ][ ][ ] [ −t ]
e−tk 0 e 0 2 3 eτ 0 e k 0
= −2t dτ ·
0 e k tk−1 0 e 2τ 3 5 0 e 2τ
0 e−2tk
[ −t ] ∫ tk [ 2τ ] [ ]
e k 0 2e 3e3τ e−tk 0
= −2t dτ ·
0 e k tk−1 3e 3τ 5e 4τ 0 e−2tk
[ −t ] [ 2τ ]t [ −t ]
e k 0 e e3τ k e k 0
=
0 e−2tk e3τ 45 e4τ t 0 e−2tk
[ ]
k−1
Finalmente se llega a
[ ]
1 − e−2T 1 − e−3T
Qk−1 = (4.32)
1 − e−3T 4 (1 − e
5 −4T )
P = Fk−1 PF⊤
k−1 + Qk−1 (4.33)
Efectuando la multiplicación e igualando los elementos individuales de las matrices
en el lado derecho e izquierdo de esta ecuación da
P11 = e−2T P11 + 2T
P12 = e−3T P12 + 3T
P22 = e−4T P22 + 5T
Resolviendo estas ecuaciones se llega a
2T
P11 =
1 − e−2T
3T
P12 =
1 − e−3T
5T
P22 =
1 − e−4T
El proceso y la medida tienen ruidos {vk } y {wk } los cuales son blancos,
con media cero, no correlacionados y además, tienen matrices de covarianza
conocidas Qk y Rk , respectivamente:
vk ∼ (0, Qk )
wk ∼ (0, Rk )
E {vk vj⊤ } = Qk δk−j (4.35)
E {wk wj⊤ } = Rk δk−j
E {vk wj⊤ } = E {wk vj⊤ } = 0
Si se tienen todas las medidas anteriores pero sin incluir el instante k disponi-
ble, para uso en la obtención del estimado de xk , entonces se puede conformar
4.4. Deducción del filtro de Kalman discreto 119
x̂−
k = E {xk |y1 , y2 , . . . , yk−1 } = estimado a priori (4.37)
x̂−
k = estimado de xk antes que se procese la medida en el instante k,
+ (4.38)
x̂k = estimado de xk después que se procese la medida en el instante k.
donde n es algún entero positivo cuyo valor depende del problema especı́fico
que se desea resolver. Si se quiere encontrar la mejor predicción de xk con
más de un paso de tiempo por delante de las medidas disponibles, entonces
se puede formar una estimación predictiva. Una manera de formar el estado
estimado predictivo es calcular el valor esperado de xk , condicionado a todas
las mediciones que estén disponibles, es decir,
donde m es algún entero positivo cuyo valor depende del problema especı́fico
que se desea resolver. La relación entre los estimados de los estados suavizado,
a posteriori, a priori y predictivo se muestra en la Figura 4.2.
En la notación que sigue, se utiliza x̂+0 para denotar el estimado inicial de
x0 antes de que sea disponible cualquier medida. La primera medida se toma
en el instante k = 1. Puesto que no se tiene ninguna medida disponible para
120 Capı́tulo 4. Filtro de Kalman discreto
Estimado suavizado
Estimado a posteriori
Estimado a priori
Predicción
1 2 3 4 5 6 7 8 9 tiempo
Figura 4.2.: Lı́nea temporal que muestra la relación entre los estimados de estado
suavizado, a posteriori, a priori y predictivo. En esta figura, se supo-
ne que se han recibido las medidas, en pasos de tiempo hasta k = 5,
inclusive. Una estimación del estado en k < 5 se denomina una estima-
ción suavizada. Una estimación del estado en k = 5 se denomina una
estimación a posteriori. Una estimación del estado en k = 6 se llama
la estimación a priori. Una estimación del estado en k > 6 se llama
predicción.
0 = E {x0 }
x̂+ (4.41)
P− − − ⊤
k = E {(xk − x̂k )(xk − x̂k ) },
+ ⊤
(4.42)
k = E {(xk − x̂k )(xk − x̂k ) }.
P+ +
Estas relaciones se explican en la Figura 4.3. La figura muestra que una vez
hecha la medida en el instante (k − 1) se obtiene un estimado de xk−1 (defini-
do como x̂+ k−1 ) y la covarianza del error de estimación (definida como Pk−1 ).
+
Cuando llega el instante k, antes que se procesen las medidas en ese instante,
se calcula un estimado de xk (denotado como x̂− k ) y la covarianza del error de
estimación (denotada como P− k ). Entonces se procesa la medida en el instante
k para refinar el estimado de xk . El valor estimado resultante de xk se deno-
ta como x̂+ +
k y su covarianza como Pk . Se comienza el proceso de estimación
con x̂+ +
0 , el mejor estimado del estado inicial x0 . Dado x̂0 , ¿cómo se podrá
− −
calcular x̂1 ? Para ello se establece x̂1 = E {x1 }. Pero nótese que x̂+
0 = E {x0 },
4.4. Deducción del filtro de Kalman discreto 121
Figura 4.3.: Lı́nea de tiempo que muestra los valores estimados de los estados
y la covarianza del error de estimación a priori y a posteriori,
respectivamente.
x̂− +
1 = F0 x̂0 + G0 u0 (4.43)
x̂− +
k = Fk−1 x̂k−1 + Gk−1 uk−1 (4.44)
en el estimado inicial de x0 :
⊤
0 = E {(x0 − x̄0 )(x0 − x̄0 ) }
P+
+ ⊤
(4.45)
= E {(x0 − x̂+
0 )(x0 − x̂0 ) }
−
Dado el valor inicial de la covarianza, P+
0 , ¿cómo calcular el siguiente valor, P1 ?
De la Ecuación (4.4) se puede recordar la forma como se propaga la covarianza
del estado de un sistema lineal discreto en el tiempo: Pk = Fk−1 Pk−1 F⊤ k−1 +
Qk−1 . De aquı́, se obtiene:
P− + ⊤
1 = F0 P0 F0 + Q0 (4.46)
P− + ⊤
k = Fk−1 Pk−1 Fk−1 + Qk−1 (4.47)
Kk = Pk−1 H⊤ ⊤
k (Hk Pk−1 Hk + Rk )
−1
= Pk H⊤ −1
k Rk
x̂k = x̂k−1 + Kk (yk − Hk x̂k−1 )
(4.48)
Pk = (I − Kk Hk )Pk−1 (I − Kk Hk )⊤ + Kk Rk K⊤
k
= (P−1 ⊤ −1
k−1 + Hk Rk Hk )
−1
= (I − Kk Hk )Pk−1
donde x̂k−1 y Pk−1 son el valor estimado y la covarianza del error de estimación
antes que la medida yk sea procesada y x̂k y Pk son el valor estimado y la
covarianza del error de estimación después que la medida yk ha sido procesada.
En esta sección, x̂− −
k y Pk son el valor estimado y la covarianza del error de
4.4. Deducción del filtro de Kalman discreto 123
Tabla 4.1.: Relaciones entre las estimaciones y las covarianzas de las secciones
3.4 y 4.4.
Estimación por mı́nimos cuadrados Filtro de Kalman
x̂k−1 = estimación antes de procesar yk ⇒ x̂−
k = estimado a prior i
Pk−1 = covarianza antes de procesar yk ⇒ P−k = covarianza a priori
x̂k = estimación después de procesar yk ⇒ +
x̂k = estimado a posterior i
Pk = covarianza después de procesar yk ⇒ P+k = covarianza a posteriori
Kk = P− ⊤ − ⊤
k Hk (Hk Pk Hk + Rk )
−1
= P+ ⊤ −1
k Hk Rk
− −
k = x̂k + Kk (yk − Hk x̂k )
x̂+
− ⊤ ⊤ (4.49)
k = (I − Kk Hk )Pk (I − Kk Hk ) + Kk Rk Kk
P+
= [(P−
k)
−1
+ H⊤ −1
k Rk Hk ]
−1
= (I − Kk Hk )P−
k
0 = E {x0 }
x̂+
+ ⊤
(4.52)
0 = E {(x0 − x̂0 )(x0 − x̂0 ) }
P+ +
Para k ∈ {1, 2, · · · , ∞}, las ecuaciones de predicción del filtro de Kalman son:
x̂− +
k = Fk−1 x̂k−1 + Gk−1 uk−1 = estado estimado a priori
(4.53)
P− + ⊤
k = Fk−1 Pk−1 Fk−1 + Qk−1
⊤ −1
= P+k Hk Rk
− −
k = x̂k + Kk (yk − Hk x̂k ) = estado estimado a posteriori
x̂+
− ⊤ ⊤
(4.54)
k = (I − Kk Hk )Pk (I − Kk Hk ) + Kk Rk Kk
P+
= [(P−k)
−1
+ H⊤ −1
k Rk Hk ]
−1
= (I − Kk Hk )P−k
Observaciones:
La primera expresión para P+k anterior, se denomina la versión estabiliza-
da de Joseph de la ecuación de corrección de la medida de la covarianza.
Fue formulada en los años sesenta. Se puede demostrar que es más estable
y robusta que la tercera. Esta primera expresión garantiza que P+ k será
siempre simétrica definida positiva si P−k es simétrica definida positiva.
+
La tercera expresión para Pk es computacionalmente más simple que
4.4. Deducción del filtro de Kalman discreto 125
min E {x̃⊤
k Sk x̃k } (4.56)
x̂− +
k+1 = Fk x̂k + Gk uk (4.58)
x̂− − −
k+1 = Fk [x̂k + Kk (yk − Hk x̂k )] + Gk uk
(4.59)
= Fk (I − Kk Hk )x̂−
k + Fk Kk yk + Gk uk
P− + ⊤
k+1 = Fk Pk Fk + Qk (4.60)
P− − − ⊤
k+1 = Fk (Pk − Kk Hk Pk )Fk + Qk
= Fk P− ⊤ − ⊤
k Fk − Fk Kk Hk Pk Fk + Qk
= Fk P− ⊤ − ⊤ − ⊤ −1 − ⊤
k Fk − Fk Pk Hk (Hk Pk Hk + Rk ) Hk Pk Fk + Qk (4.61)
⊤
(4.62)
k = (I − Kk Hk )(Fk−1 Pk−1 Fk−1 + Qk−1 )
P+ +
[ ] [ ][ ]
ṙ 0 1 0 r
v̇ = 0 0 1 v
ȧ 0 0 0 a
ẋ = Ax
x̂− +
k = Fx̂k−1
P− + ⊤
k = FPk−1 F + Qk−1 con Qk−1 = 0
⊤
= FP+
k−1 F
Kk = P− ⊤ − ⊤
k Hk (Hk Pk Hk + Rk )
−1
Si se escribe la matrix P− k ∈ R
3×3 , en términos de sus elementos, dado que es
simétrica, y se sustituye para Hk y Rk en la ecuación anterior, se obtiene
− − [ ] − − [ ]
−1
p11,k p−
12,k p13,k 1 p11,k p−
12,k p13,k 1
Kk = p− − −
21,k p22,k p23,k 0 [1 0 0]p− − −
21,k p22,k p23,k 0 + σ2
p− − −
31,k p32,k p33,k
0 p− − −
31,k p32,k p33,k
0
− −
p11,k ( )−1 ( )−1 p11,k
Kk = p−
21,k
p− + σ 2
11,k ≡ p−
11,k + σ
2 p−
12,k
− −
p31,k p13,k
−
p−
11,k p11,k 0 0
− 1 p− − − 1 p− −
k = Pk − −
P+ 2 12,k [1 0 0] Pk = Pk − − 2 12,k 0 0 Pk
p11,k + σ p11,k + σ
p−
13,k p−
13,k 0 0
− − − − 2 −
p11,k p12,k p13,k (p11,k ) − −
p11,k p12,k p− 11,k p13,k
1
= p− − −
12,k p22,k p23,k − − 2
p− −
11,k p12,k (p−12,k )
2 p− −
12,k p13,k
p + σ
p− − −
13,k p23,k p33,k
11,k p− − − −
11,k p13,k p12,k p13,k (p−
13,k )
2
Se usará esta expresión para mostrar que del tiempo k − al tiempo k + la traza
de la covarianza del error de estimación decrece.
132 Capı́tulo 4. Filtro de Kalman discreto
3500
P−2
1500
1000
P+ P+ P5
+
3 4
500
P+2
P+ +
P1
0 0
0 1 2 3 4 5
Pasos de tiempo
Figura 4.5.: Los primeros cinco pasos de las varianzas en la estimación del
error a priori y a posteriori del Ejemplo 4.5.
Puesto que la matriz de covarianza es simétrica, entonces, de la Ecuación (??)
se puede ver que la traza de P+k , está dada como
( ) ( ) ( )
(p−
11,k )
2 (p−
12,k )
2 (p−
13,k )
2
Tr(P+
k) = p− − + p− − + p− −
11,k
p−
11,k + σ
2 22,k
p−
11,k + σ
2 33,k
p−
11,k + σ
2
(p− 2 − 2 −
11,k ) + (p12,k ) + (p13,k )
2
= Tr(P−
k) − (4.63)
p−
11,k + σ
2
Cuando se toma una nueva medida, se espera que el estimado del estado mejore.
Es decir, se espera que la covarianza disminuya, y la ecuación anterior muestra
que ésta en realidad ha decrecido. Esto es, la traza de P+ k es menor que la tra-
−
za de Pk . Este sistema fue simulado con cinco unidades de tiempo entre pasos
de discretización (T = 5), y una desviación estándar posición–medida de 30 uni-
dades. La Figura 4.5 muestra la varianza de la estimación de la posición (p− 11,k
y p+
11,k ) para los primeros cinco pasos de tiempo del filtro Kalman. Se puede ob-
servar que la varianza (incertidumbre) se incrementa desde el estado de tiempo
uno al siguiente, pero luego decrece en cada paso de tiempo cuando la medida es
procesada.
4.5. Propiedades del filtro de Kalman 133
3500
2500
2000
1500
1000
500
0
0 10 20 30 40 50 60
Pasos de tiempo
Figura 4.6.: Los primeros sesenta pasos de las varianzas del error de estima-
ción a priori y a posteriori del Ejemplo 4.5.
Error de medición
Error de estimation
10 20 30 40 50 60
Pasos de tiempo
varianza aumenta entre pasos de tiempo, y luego decrece en cada paso de tiempo.
También se puede observar de esta figura que la varianza converge al valor de
estado estacionario.
La Figura 4.7, muestra el error posición de la medida (con una desviación
estándar de 30) y el error de estimación de la posición a posteriori. El error de la
estimación comienza con una desviación estándar cercana a 30, pero al final de la
simulación la desviación estándar está cerca de 11.
P− + ⊤
k = Fk−1 Pk−1 Fk−1 + Qk−1 (4.64)
P+
k = P−
k − P− ⊤ − ⊤
k Hk (Hk Pk Hk
−1
+ Rk ) H k P−
k
Si la matriz P−
k ∈ R
n×n
se factoriza como
P− −1
k = Ak Bk (4.65)
P− −1
k+1 = Ak+1 Bk+1 (4.66)
B−1 −⊤ ⊤ −1 −⊤
k+1 = [Fk Hk Rk Hk Ak + Fk Bk ]
−1
= [F−⊤ ⊤ −1 −1
k (Hk Rk Hk Ak Bk + I)Bk ]
−1
= B−1 ⊤ −1 −1 −1 ⊤
k [Hk Rk Hk Ak Bk + I] Fk (4.69)
Ak+1 B−1 −⊤ ⊤ −1 −⊤
k+1 = [(Fk + Qk Fk Hk Rk Hk )Ak + Qk Fk Bk ]×
× B−1 ⊤ −1 −1 −1 ⊤
k [Hk Rk Hk Ak Bk + I] Fk
Ak+1 B−1 −⊤ ⊤ −1 −1 −⊤
k+1 = [(Fk + Qk Fk Hk Rk Hk )Ak Bk + Qk Fk ]
× [H⊤ −1 −1 −1 ⊤
k Rk Hk Ak Bk + I] Fk (4.70)
Sustituyendo Ak B−1 −
k por Pk en la ecuación anterior da
Ak+1 B−1 −⊤ ⊤ −1 − −⊤ ⊤ −1 − −1 ⊤
k+1 = [(Fk + Qk Fk Hk Rk Hk )Pk + Qk Fk ][Hk Rk Hk Pk + I] Fk
= [Fk P− −⊤ ⊤ −1 − −⊤ ⊤ −1 − −1 ⊤
k + Qk Fk Hk Rk Hk Pk + Qk Fk ][Hk Rk Hk Pk + I] Fk
= [Fk P− −⊤ ⊤ −1 − ⊤ −1 − −1 ⊤
k + Qk Fk (Hk Rk Hk Pk + I)][Hk Rk Hk Pk + I] Fk
= F k P− ⊤ −1 − −1 ⊤
k [Hk Rk Hk Pk + I] Fk + Qk (4.71)
Ak+1 B−1 − ⊤ − ⊤ −1 − ⊤
k+1 = Fk Pk [I − Hk (Hk Pk Hk + Rk ) Hk Pk ]Fk + Qk
= Fk [P− − ⊤ − ⊤ −1 − ⊤
k − Pk Hk (Hk Pk Hk + Rk ) Hk Pk ]Fk + Qk
⊤
= Fk P+
k Fk + Qk
= P−
k+1 (4.72)
136 Capı́tulo 4. Filtro de Kalman discreto
K∞ = P − ⊤ − ⊤ −1
∞ H (HP∞ H + R) , la ganancia del filtro de Kalman en estado
estacionario.
Sistemas escalares
Se puede usar la Ecuación (4.68) para obtener una solución en forma cerrada
del filtro escalar de Kalman para sistemas invariantes en el tiempo. Supóngase
que F, Q, H y R son escalares constantes. Entonces de la Ecuación (4.68) se
obtiene
[ ] [ 2
][ ] [ ]
Ak+1 F + HF RQ Q
F
Ak Ak
= H 2 1
=Ψ (4.78)
Bk+1 FR F
Bk Bk
Nótese la secuencia
[ ] [ ] [ 2 ]
λ1 0 −1 λ1 0 −1 λ1 0
2
Ψ =M M ×M M =M M−1
0 λ2 0 λ2 0 λ22
[ 2 ] [ ] [ 3 ]
λ1 0 −1 λ1 0 −1 λ1 0
3
Ψ =M M ×M M =M M−1
0 λ22 0 λ2 0 λ32
[ ]
λk1 0
k
Ψ =M M−1 (4.80)
0 λk2
Ejemplo 4.6
Propagación de una covarianza escalar. Considérese el siguiente sistema escalar:
xk+1 = xk + vk
y k = x k + wk
vk ∼ (0, 1), wk ∼ (0, 1)
Este es un ejemplo muy simple pero que puede surgir en muchas aplicaciones;
v.gr.: puede representar algún parámetro xk con variación lenta que se mide direc-
tamente. El término de ruido del proceso vk se tiene en cuenta en las variaciones
de xk y, la medida del término de ruido wk , se considera para errores en la medida.
En este sistema, se tiene F = H = Q = R = 1.
Sustituyendo estos valores en la Ecuación (4.82) da
√
µ1 = 3 + 5
√
µ2 = 3 − 5
√
ν1 = 1 + 5
√
ν2 = 1 − 5
− −
1 (2P1 − ν2 ) − ν2 µ2 (2P1 − ν1 )
ν1 µk−1 k−1
Pk− = − −
1 (2P1 − ν2 ) − 2µ2 (2P1 − ν1 )
2µk−1 k−1
2
Ganancia de Kalman
Covarianza del error de estimación
1.8
1.6
1.4
1.2
0.8
0.6
0 1 2 3 4 5 6 7 8 9 10
tiempo
La teorı́a presentada en este capı́tulo hace del filtro de Kalman una escogencia
atractiva para la estimación de estados. Pero cuando un filtro de Kalman se
implementa en un sistema real éste puede no funcionar, aun cuando la teorı́a
sea correcta. Dos de las primeras causas para la falla del filtrado de Kalman
son: la aritmética de precisión finita y los errores de modelado. Esta metodo-
logı́a supone que la aritmética del filtro de Kalman tiene precisión infinita. En
4.5. Propiedades del filtro de Kalman 141
Ejemplo 4.7
A continuación se va a ilustrar el uso de ruido de un proceso ficticio. Supóngase que
se está tratando de estimar un estado que se asume constante, pero en realidad es
una rampa. En otras palabras, se tiene un error de modelado. El modelo asumido
(pero incorrecto), sobre el cual se basa el filtro de Kalman, se da de la siguiente
manera:
xk+1 = xk + vk
yk = xk + wk
(4.84)
vk ∼ (0, 0)
wk ∼ (0, 1)
El ruido del proceso se asume cero, lo que significa que se está modelando xk
como una constante. De la Ecuación (4.54) se obtienen las ecuaciones del filtro de
142 Capı́tulo 4. Filtro de Kalman discreto
400
300
200
100
0
0 10 20 30 40 50
Tiempo
600
Estado verdadero
x estimado (Q = 1)
500 x estimado (Q = 0.1)
x estimado (Q = 0.01)
x estimado (Q = 0)
400
300
200
100
0
0 10 20 30 40 50
Tiempo
0.7
0.6
Q = 1
Q = 0.1
0.5
Q = 0.01
Q = 0
0.4
0.3
0.2
0.1
0
0 10 20 30 40 50
Tiempo
Figura 4.11.: Ganancia de Kalman para varios valores del ruido del proceso.
A medida que el ruido del proceso ficticio se hace más grande, el error de estima-
144 Capı́tulo 4. Filtro de Kalman discreto
ción se hace más pequeño. Por supuesto, esto es al precio de empeorar el desem-
peño en el caso de que el modelo del sistema asumido sea realmente correcto.
El diseñador tiene que añadir una cantidad apropiada de ruido del proceso fic-
ticio para compensar el desempeño bajo condiciones nominales de mal modelado.
La Figura 4.11 muestra la evolución temporal de la ganancia de Kalman Kk
en este ejemplo, para varios valores de Q. Como era de esperar, la ganancia Kk
converge a un valor de estado estacionario más grande cuando Q es más grande,
haciendo que el filtro sea más sensible a las mediciones [ver la expresión x̂+k en
la Ecuación (4.85)]. Esto compensa los errores de modelado. También téngase en
cuenta en la Figura 4.11, que la ganancia de Kalman en estado estacionario es de
aproximadamente 0.62 cuando Q = 1.
Este ejemplo ilustra el principio general de que el ruido del modelo es bueno,
pero sólo hasta cierto punto. Si un modelo de sistema tiene mucho ruido, enton-
ces es difı́cil estimar sus estados. Pero si un modelo de sistema tiene muy poco
ruido, entonces el estimador de estados puede ser demasiado susceptible a errores
de modelado. En el diseño de un modelo para un filtro de Kalman, se necesita
balancear la confianza en el modelo (bajo nivel de ruido resulta en un seguimiento
cercano del modelo, es decir, bajo ancho de banda) con una autoduda saludable
(alto ruido resulta en mayor capacidad de respuesta del filtro, es decir, alto ancho
de banda).
Un examen a las ecuaciones del filtro muestra por qué la adición de ruido al
proceso ficticio compensa los errores de modelado. Recordando las ecuaciones
del filtro de Kalman, (4.53) y (4.54), algunas de las cuales se repiten aquı́:
4.6. Resumen
En este capı́tulo, se ha presentado la esencia del filtro Kalman en tiempo dis-
creto. Durante las últimas décadas, este algoritmo de estimación ha encontrado
aplicaciones en prácticamente todas las áreas de la ingenierı́a.
Se ha mostrado que las ecuaciones del filtro de Kalman se pueden escribir de
varias maneras diferentes, cada una de las cuales puede parecer muy diferente
a las otras, aunque todas ellas son matemáticamente equivalentes. Se ha visto
que el filtro de Kalman es óptimo incluso cuando el ruido no es gaussiano. El
filtro de Kalman es el estimador óptimo cuando el ruido es gaussiano, y es el
estimador lineal óptimo cuando el ruido no es gaussiano.
También se ha observado que el filtro de Kalman puede no funcionar adecua-
damente si las suposiciones subyacentes no se mantienen. Para esta situación
se mencionaron brevemente algunas formas de compensar los supuestos viola-
dos. En los capı́tulos siguientes se expandirán y generalizarán los resultados
presentados en éste.
Bibliografı́a
Gibbs, B. P. (2011). Advanced Kalman filtering, least-squares and modeling. John Wiley &
Sons, Inc., Hoboken, New Jersey, EE UU.
Grewal, M. S. and Andrews, A. P. (2008). Kalman Filtering: Theory and Practice Using
MATLAB. John Wiley & Sons, Inc., Hoboken, New Jersey, EE UU., 3a edition.
Kailath, T. (1968). An innovations approach to least-squares estimation: Part 1. linear
filtering in additive white noise. IEEE Trans. Autom. Control, AC-13(6): 646–655.
Kailath, T., Sayed, A., and Hassibi, B. (2000). Linear Estimation. Upper Saddle River, New
Jersey, USA.
Kwakernaak, H. and Sivan, R. (1972). Linear Optimal Control Systems. John Wiley & Sons,
Inc., New York, USA.
Åström, K. J. and Wittenmark, B. (1997). Computer-Controlled Systems: Theory and De-
sign. Prentice Hall Information and System Sciences, 3a edition.
147
Capı́tulo 5
5.1. Introducción
149
150 Capı́tulo 5. Formas alternas del filtro de Kalman
Número de condición
Recordando la definición de valores singulares de una matriz. Una matriz P
n × n tiene n valores singulares σ dados como
σ 2 (P) = λ(P⊤ P)
= λ(PP⊤ ) (5.1)
Algunos autores usan definiciones alternas para número de condición, por ejem-
plo, pueden definir el número de condición de una matriz como el cuadrado de
la anterior definición. Cuando κ(P) → ∞, la matriz P se dice ser pobremente
condicionada o mal condicionada, y P se aproxima a una matriz singular. En
la implementación de un filtro de Kalman, la matriz
{ de covarianza} del error
⊤
P será siempre definida positiva porque P = E (x − x̂)(x − x̂) . Se usa la
notación estándar
P>0 (5.3)
para indicar que P es definida positiva. Esto implica que P es invertible, lo cual
a su vez es equivalente a decir que todos los autovalores de P son mayores que
cero. Pero supóngase en el filtro de Kalman que algunos de los elementos de x
son estimados a mucha más precisión que otros elementos de x. Por ejemplo,
supóngase que
[ 6 ]
10 0
(5.4)
0 10−6
for j = 1 to n do
γji = 0 j<i
1 ( ∑ )
γji = pji − i−1 γ γ
k=1 jk ik j>i
γii
end for
end for
Ejemplo 5.1
Supóngase que se tiene una matriz P dada como
[ ]
1 2 3
P= 2 8 2
3 2 14
El algoritmo de factorización de Cholesky dice que, para i = 1,
√
γ11 = p11 = 1
1
γ21 = (p21 ) = 2
γ11
1
γ31 = (p31 ) = 3
γ11
Para i = 2, el algoritmo dice que
γ12 = 0
√ √
∑1 √
γ22 = p22 − 2 =
γ2j p22 − γ21
2 = 8−4=2
j=1
( )
1 ∑1
1 1
γ32 = p32 − γ3k γ2k = (p32 − γ31 γ21 ) = (2 − 3 × 2) = −2
γ22 γ22 2
k=1
Ası́ se obtiene
[ ]
1 0 0
Γ= 2 2 0
3 −2 1
Recordando que para una matriz general A se tiene λ(A2 ). Por lo tanto, se ve
desde la anterior ecuación que
σ 2 (P) = [σ 2 (Γ)]2
σmax (P) σ 2 (Γ)
= max2
(5.8)
σmin (P) σmin (Γ)
2
κ(P) = κ (Γ)
θ1⊤ θ2 = θ2⊤ θ1 = 0
(5.14)
θ1⊤ θ1 = θ2⊤ θ2 = I
Se puede usar esta ecuación, junto con la ecuación (5.14), para escribir
⊤/2
Γ− − ⊤ ⊤ ⊤ ⊤ ⊤
+ + 1/2
k (Γk ) = Fk−1 Γk−1 θ1 θ1 (Γk−1 ) Fk−1 + Qk−1 θ2 θ2 Qk−1
⊤ ⊤ 1/2 ⊤/2
= Fk−1 Γ+ +
k−1 (Γk−1 ) Fk−1 + Qk−1 Qk−1 (5.15)
156 Capı́tulo 5. Formas alternas del filtro de Kalman
Si Γ+ +
k−1 es la raı́z cuadrada de Pk−1 , esto implica que
P− + ⊤
k = Fk−1 Pk−1 Fk−1 + Qk−1 (5.16)
El algoritmo de Householder
(k)
donde Aik es el elemento de la i−ésima fila y la k−ésima columna
de A(k) . La función sgn(x) se define como
{
+1 x≥0
sgn =
−1 x<0
3. Después de que los anteriores pasos han sido ejecutados, A(n+1) tiene la
forma
[ ]
(n+1) W
A = (5.24)
0
1
βk = ( )
(k)
ςk ςk + Akk
for i = 1 to 2n do
if i < k then
(k)
ui = 0
else if i = k then
(k) (k)
ui = ςk + Akk
else if i > k then
(k) (k)
ui = Aik
end if
end for
for i = 1 to n do
if i < k then
(k)
yi = 0
else if i = k then
(k)
yi = 1
else if i > k then
(k) (k)
yi = βk u(k)⊤ Ai
end if
end for
[ ]
W
A (k+1)
=A (k)
−u y (k) (k)⊤
→A (n+1)
=
0
for i = 1 to n do
∏1I −(i)βk u u
θ (k) = (k) (k)⊤
θ = nθ
end for
end for
160 Capı́tulo 5. Formas alternas del filtro de Kalman
Ejemplo 5.2
Usar el método de Householder para encontrar una
matriz
ortogonal θ tal que
[ ] 1 1
W 2 2
donde W es una matriz 2 × 2 y A = = A(1)
1
θA =
0 0
2 2
Solución:
Aplicando el Algoritmo 5.2 se obtiene:
v
( )u
u∑4 ( ) √
(1) t (1) 2
ς1 = sgn A11 Ai1 = 12 + 2 2 + 0 2 + 2 2 = 3
i=1
1 1 1
β1 = ( )= =
ς1 ς1 +
(1)
A11 3(3 + 1) 12
(1) (1)
u1 = ς1 + A11 = 3 + 1 = 4
(1) (1)
u2 = A21 = 2
(1) (1)
u3 = A31 = 0
(1) (1)
u4 = A41 = 2
[ ]⊤
⊤
u(1) = u(1)
1
(1)
u2
(1)
u3
(1)
u4 = [4 2 0 2]
(1)
y1 = 1
1
(1) (1) 1 2 1
y2 = β1 u(1)⊤ A2 = [4 2 0 2] = (4 + 4 + 0 + 4) = 1
12 1 12
2
⊤
y = [1 1]
1 1 4 1 1 4 4
2 2 2 2 2 2 2
A(2) = A(1) − u(1) y(1)⊤ = − −
1 0
[1 1] =
0 0 1 0 0
2 2 2 2 2 2 2
−3 −3
0 0
A(2) =
0 1
0 0
5.2. Filtrado de raı́z cuadrada 161
v
u 4 ( ) √( ) ( ) ( )
u∑ (2) 2 (2) 2 (2) 2 (2) 2
ς2 = sgn(0) t Ai2 = A22 + A32 + A42
i=2
√
= 0+1+0=1
1 1
β2 = ( )= =1
(2)
ς2 ς2 + A22 1(1 + 0)
(2)
u1 = 0
(2)
u2 = ς2 + 0 = 1
(2) (2)
u3 = ς2 + A32 = 1 + 0 = 1
(2) (2)
u4 = A42 = 0
0
(2) 1 ⊤
u = = [0 1 1 0]
1
0
(2)
y1 = 0
(2)
y1 = 0
⊤
y(2) = [0 1]
−3 −3 0 −3 −3 0 0
0 0 1 0 0 0 1
A(3) = A(2) − u(2) y(2)⊤ = − −
1 1
[0 1] =
0 0 1 0 1
0 0 0 0 0 0 0
−3 −3
0 −1
A(3) =
0 0
0 0
[ ]
−3 −3
W= (5.26)
0 −1
162 Capı́tulo 5. Formas alternas del filtro de Kalman
( )( )
θ = I − β2 u(2) u(2)⊤ I − β1 u(1) u(1)⊤
1 0 0 0 0 1 0 0 0 4
0 1 0 0 1 0 1 0 0 1 2
= −1 −
0 0 1 0 1
[0 1 1 0]
0 120
[4 2 0 2]
0 0 1
0 0 0 1 0 0 0 0 1 2
1 1
1 0 0 0 −3 − 23 0 − 32 −3 − 23 0 − 32
0 0 −1 0 − 3
2 2
0 − 13 0
= 0 −1 0
θ=
3
0 −1 0 0 0 0 1 0 23 − 23 0 1
3
0 0 0 1 −2 3 − 13 0 32 − 23 − 13 0 2
3
(k)
donde Ak es la k-ésima columna de A(k) .
5.2. Filtrado de raı́z cuadrada 163
∑
k−1
( )
θk := θk − θk θi⊤ θi
i=1
θk
θk := (5.33)
||θ||2
v
u
√ u 1
u 2 √
=u
(1)⊤ (1)
ς1 = A1 A1 t[1 2 0 2] 0 = 1 + 4 + 0 + 4 = 3
2
1
1 (1)⊤ (1) 1 2
W11 = 3; W12 = A1 A1 = [1 2 0 2] = 3
3 3 0
2
W1 = [W11 W12 ] = [3 3]
1 (1)⊤ 1 [ ]
θ1 = A 1 = [1 2 0 2] = 13 2
3 0 2
3
ς 3
5.2. Filtrado de raı́z cuadrada 165
1 1 0
(2) (1) 1 (1) 2 1 2 0
A2 = A2 − W12 A1 = − ·3· =
ς1 1 3 0 1
2 2 0
v
u
√ u 0
u 0 √
ς2 =
(2)⊤
A2
(2)
A2 u
= t[0 0 1 0] = 0 + 0 + 1 + 0 = 1
1
0
Para k = 2
1 (2)⊤ 1
θ2 = A = [0 0 1 0] = [0 0 1 0]
ς2 2 1
θ3 = [1 0 0 0]
θ4 = [0 1 0 0]
θ5 = [0 0 1 0]
θ6 = [0 0 0 1]
( ) ( )
θ3 = θ3 − θ3 θ1⊤ θ1 − θ3 θ2⊤ θ2
1
3
2 [ 1 ]
= [1 0 0 0] −
[1 0 0 0] 3
0 3
2
3 0 2
3 −0
2
3
1 [1 2 ]
= [1 0 0 0] − 3 3 0 23
[ 3]
θ3 = 89 − 29 0 − 92
θ3 [ √ ]
θ3 = = 2 3 2 − 3√ 1
0 − 1
√
||θ3 ||2 2 3 2
[ √ √ √ ]
θ3 = 2 3 2 − 62 0 − 62
5.2. Filtrado de raı́z cuadrada 167
( ) ( ) ( )
θ4 = θ4 − θ4 θ1⊤ θ1 − θ4 θ2⊤ θ2 − θ4 θ3⊤ θ3
1
3 0
2 [1 ] 0
= [0 1 0 0]−[0 1 0 0] 3
· 2
0 3 −[0
2 1 0 0] ·[0 0 1 0]
0 3 3 1
2 0
3
√
2 2
3√2 [ √ √ √ ]
− [0 1 0 0] − 6 · −2 2
0 2
− 62
0√ 3 6
− 62
[ ] [ ] [ ]
= [0 1 0 0] − 29 49 0 49 + 92 − 18 1
0 − 18
1
= 0 1
2 0 − 12
θ4
θ4 =
||θ4 ||2
√ √ √
1 1 1 1 2
||θ4 ||2 = 0+ +0+ = =√ =
4 4 2 2 2
[ √ √ ]
θ4 = 0 22 0 − 22
1 2 2
3 3 0 3
√ 0 1√ 0 0√
θ = 2 2
−√ 62 0 − √62
3
0 2
2
0 − 22
168 Capı́tulo 5. Formas alternas del filtro de Kalman
Ejemplo 5.4
Supóngase que en el tiempo (k − 1) el filtro de Kalman tiene la matriz del sis-
tema, la covarianza del ruido del proceso y la raı́z cuadrada de la covarianza de
estimación a posteriori dadas por
[ ] [ ] [ ]
1 1 0 0 + 1 0
Fk−1 = ; Qk−1 = ; Γk−1 =
0 1 0 2 0 1
1/2 ⊤/2
Puede verificarse que la raı́z cuadrada de Qk−1 (tal que Qk−1 Qk−1 = Qk−1 ) está
dada por
[ ]
1/2 0 0
Qk−1 =
−1 −1
⊤
Γ+ +⊤
i−1,k Γi−1,k Hik
Kik = ⊤
(5.36)
Hik Γ+ +⊤
i−1,k Γi−1,k Hik + Rik
y P+
ik puede escribirse como
( )
+ +⊤ ⊤
Γ Γ H
i−1,k i−1,k ik H ik
P+
ik = I− ⊤
Γ+ +⊤
i−1,k Γi−1,k
Hik Γ+ Γ +⊤
H
i−1,k i−1,k ik + R ik
( )
+ +⊤ ⊤ +
Γ Γ H
i−1,k i−1,k ik H ik Γ i−1,k
= Γ+ i−1,k − ⊤
Γ+⊤
i−1,k
Hik Γ+ Γ+⊤
H
i−1,k i−1,k ik + R ik
⊤
ik = Γi,k Γi,k = Γi−1,k (I − αϕϕ )Γi−1,k
P+ + +⊤ + +⊤
(5.37)
170 Capı́tulo 5. Formas alternas del filtro de Kalman
⊤ 1
ϕ = Γ+⊤
i−1,k Hik α = (5.38)
ϕ⊤ ϕ + Rik
De la Ecuación (5.38)
αϕ⊤ ϕ + αRik = 1
αϕ⊤ ϕ = 1 − αRik (5.45)
Sustituyendo (5.45) en (5.44) se obtiene
0 = (1 − αRik )γ 2 − 2γ + 1; es decir,
αRik γ 2 = 1 − 2γ + γ 2 = (1 − γ)2
√
± αRik γ = 1 − γ
√
γ ± αRik γ = 1
1
γ= √ (5.46)
1 ± αRik
El cual es el factor de ajuste de modo que se cumpla la Ecuación (5.39).
172 Capı́tulo 5. Formas alternas del filtro de Kalman
2. Realizar
for i = 1 to r do
(a) Definir:
(i) Hik ← i-ésima fila de Hk ,
(ii ) yik ← i-ésimo elemento de yk ,
(iii ) Rik ← varianza de la i -ésima medición (Rk diagonal).
(b) Calcular:
⊤
ϕi = Γ+⊤
i−1,k Hik
1
αi = ⊤
ϕi ϕi + Rik
1
γ= √ (5.48)
1 ± αi Rik
⊤
i−1,k (I − αi γi ϕi ϕi )
Γik = Γ+
+
Kik = αi Γ+ ik ϕi
x̂ik = x̂i−1,k + Kik (yik − Hik x̂+
+ +
i−1,k )
end for
3. raı́z cuadrada de la covarianza a posteriori y el estimado del estado a
posteriori
Γ+ +
k = Γrk
(5.49)
x̂+ +
k = x̂rk
5.2. Filtrado de raı́z cuadrada 173
Ejemplo 5.6
Supóngase que se tiene un sistema LTI con
[ ] [ ]
1 0 1 0
P− = F =
k 0 1 0 1
[ ]
0 0
H = [1 0] Q=
0 0
sólo se tiene que iterar la Ecuación (5.48) una vez, ya que sólo se tiene una
medición. Los valores alrededor de los parámetros dados en la Ecuación (5.48)
son
[ ]
− ⊤ ⊤ 1
ϕ = (Γk ) H =
0
1 1
α= ⊤ = ≈1
ϕ ϕ+R 1+R
1 1
γ= √ = √
1 + αR 1+ R
[ √ ]
√R
− ⊤ 0
k = Γk (I − αγϕϕ ) =
Γ+ 1+ R
0 1
Nótese que Γ+ +⊤
k Γk es no singular. Los valores alrededor de la raı́z cuadrada de la
covarianza a priori, los parámetros de la Ecuación (5.48) y la ganancia de Kalman
en el siguiente paso de tiempo (k + 1) estarán dados por
Γ− +
k+1 = Γk
[ √ ]
√R
ϕ = (Γ− ⊤ ⊤
k+1 ) H =
1+ R
0
√ √
1 1+R+2 R 1+2 R (5.53)
α= ⊤ = √ ≈ √
ϕ ϕ+R R2 + 2R + 2R R 2R(1 + R)
√ [ ] [ 1√
]
R √
− 1+2 R
Kk+1 = αΓk+1 ϕ ≈ √ 1+R+2 R ≈ 2(1+ R)
2R(1 + R) 0 0
I = P−1 (5.55)
P− + ⊤
k = Fk−1 Pk−1 Fk−1 + Qk−1 (5.58)
I− −1 −1 ⊤ −1 −1 ⊤ −1
k = Qk−1 − Qk−1 Fk−1 (Ik−1 + Fk−1 Qk−1 Fk−1 ) Fk−1 Qk−1
+
(5.61)
I− −1 −1 ⊤ −1 −1 ⊤ −1
k = Qk−1 − Qk−1 Fk−1 (Ik−1 + Fk−1 Qk−1 Fk−1 ) Fk−1 Qk−1
+
− ⊤ −1
Ik = Ik + Hk Rk Hk
+
−1 ⊤ −1
Kk = (I+ k ) Hk Rk (5.62)
− +
x̂k = Fk−1 x̂k−1 + Gk−1 uk−1
− −
k = x̂k + Kk (yk − Hk x̂k )
x̂+
end for
la medida promedio en cada paso de tiempo? Usar las ecuaciones del filtro
estándar de Kalman para calcular las matrices de covarianza a priori y a
posteriori en k = 1 y k = 2 y comprobar que son las inversas de las matrices
de información que se calcularon en la parte 1.
Solución:
I− −1 −1 2 −1 −1 −1
k = Qk−1 − Qk−1 F (Ik−1 + F Qk−1 ) Qk−1
2 +
1 F2 1
= − ·
Qk−1 Q2k−1 I+
k−1 +
F2
Qk−1
1 F2 1
I−
1 = − 2 · + F2
Q 0 Q 0 I0 +
Q0
1 1 4
=1− · =
4 1 + 14 5
− ⊤ −1 − ⊤ −1
I+
k = Ik + Hk Rk Hk = Ik + Hk I Hk
[ ][ ]
− 1 0 1
I+
k = Ik + [1 1] 0 1 1
4 14
I+
1 = +2=
5 5
1 F 2 1
I−
2 = − 2 · + F2
Q 1 Q 1 I1 +
Q 1
4 1 42 1 4 4 1 56
= − · 2· = − · =
5 4 5 14
5 + 14 · 4
5
5 25 14
5 + 1
5
75
−
I+
2 = I2 + 2
56 206
I+
2 = +2=
75 75
5.4. Resumen 179
2. Supóngase que se tienen dos mediciones independientes del escalar x, cada una
con varianza R. Entonces la varianza de la medición promedio se calcula como
{[ ]2 } {[ ]2 }
1 1 1
E (y1 + y2 ) − x =E (y1 − x) + (y2 − x)
2 2 2
{[ }
w1 w2 ]2
=E +
2 2
R R R
= + =
4 4 2
En este problema cada medición tiene una varianza unidad, por lo que la varianza
de la medición promedio es 1/2. Se pueden utilizar las ecuaciones (4.53) y (4.54)
para obtener
Pk− = F 2 Pk−1
+
+ Qk−1
P1− = F 2 P0+ + Q0
1 5
= +1=
4 4
P1+ = [(P1− )−1 + H ⊤ R−1 H]−1
( )−1
4 5
= +2 =
5 14
1 5 5 75
P2− = · + =
4 14 4 56
− −1 −1
P2+ = [(P2 ) + 2]
( )−1
56 75
= +2 =
75 206
Nótese que los cantidades Pi− y Pi+ obtenidas aquı́, son las inversas de las canti-
dades I−1 e I1 obtenidas en la parte 1.
+
5.4. Resumen
En este capı́tulo se han estudiado algunos modelos alternos del filtro de Kal-
man. Uno de estos modelos es el filtro de información.
El filtro de información es equivalente al filtro Kalman, pero se propaga a
la inversa de la covarianza. Esto puede ser computacionalmente beneficioso en
180 Capı́tulo 5. Formas alternas del filtro de Kalman
los casos en los que el número de mediciones es mucho mayor que el número de
estados. El filtrado de raı́z cuadrada aumenta de manera efectiva la precisión
del filtro de Kalman.
Aunque estos enfoques requieren esfuerzo computacional adicional, pueden
ayudar a prevenir la divergencia y la inestabilidad. Hay una excelente y com-
pleta visión del filtrado de raı́z cuadrada en [Bierman, 1977].
Existe un número de diferentes opciones en la aplicación de un filtro de
Kalman. En este capı́tulo se han tratado dos modelos
Filtrado de covarianza o filtrado de información
de filtrado estándar, filtrado de raı́z cuadrada.
Cualquiera de estas opciones se pueden hacer de forma independiente de las
demás. Por ejemplo, se pude optar por combinar el filtrado de filtrado informa-
ción con el filtrado de raı́z cuadrada [Kaminski et al., 1971a]. Una comparación
numérica de diversas formulaciones de filtro de Kalman (incluyendo el filtro
estándar, el filtro de covarianza raı́z cuadrada, el filtro de la información de
la raı́z cuadrada, y el algoritmo de Chandrasekhar) se da en [Verhaegen and
Van Dooren, 1986]. Comparaciones numéricas y computacionales de distin-
tos enfoques de filtrado de Kalman se dan en [Bierman, 1973], [Bierman and
Thornton, 1977]. Filtrado raı́z cuadrada de tiempo continuo se discute en [Morf
et al., 1978].
Bibliografı́a
Andrews, A. (1968). A square root formulation of the Kalman covariance equations. AIAA
Journal, 6(6): 1165–1166.
Battin, R. (1964). Astronautical Guidance. Electronic sciences series. McGraw-Hill, New
York.
Bellantoni, J. and Dodge, K. (1967). A square root formulation of the Kalman-Schmidt
filter. AIAA Journal, 5(7): 1309–1314.
Bierman, G. (1973). A comparison of discrete linear filtering algorithms. IEEE Transactions
on Aerospace and Electronic Systems, AES-9(1): 28–37.
Bierman, G. and Thornton, C. (1977). Numerical comparison of Kalman filter algorithms:
Orbit determination case study. Automatica, 13(1): 23–35.
Bierman, G. J. (1977). Factorization Methods for Discrete Sequential Estimation, volume
128 of Monographs and Textbooks. Academic Press,Inc, 111 Fifth Avenue, New York.
New York 10003, NY, USA.
Brown, R. and Hwang, P. (1996). Introduction to Random Signals and Applied Kalman
Filtering. New York, USA.
Chen, C.-T. (1984). Linear System Theory and Design. Saunders College Pu., Harcourt
Brace Jovanovich, Inc., Orlando, FL 32887, USA, 2a edition.
Dyer, P. and McReynolds, S. (1969). Extension of square–root filtering to include process
noise. Journal of Optimization Theory and Applications, 3(6): 444–458.
Faddeev, D. and Faddeeva, V. (1963). Computational Methods of Linear Algebra. W. H.
Freeman and Company, San Francisco, USA.
Golub, G. and Van Loan, C. (1996). Matrix Computations. The Johns Hopkins University
Press, 2715 North Charles Street, Baltimore, Maryland 21218-4319, USA, 3a edition.
Horn, R. and Johnson, C. (1985). Matrix Analysis. Cambridge University Press, New York,
USA.
Householder, A., S. (1975). The Theory of Matrices in Numerical Analysis. Reimpresión del
libro publicado originalmente en 1964. 180 Varick Street NY 10014, NY, USA.
Jordan, T. L. (1968). Experiments on error growth associated with some linear least squares
procedures. Mathematics of Computation, AMS, (22):579–588.
Kailath, T. (1980). Linear Systems. Prentice–Hall, Inc., Englewoods Cliffs, N. J., 07632,
181
182 BIBLIOGRAFı́A
USA.
Kaminski, P., Bryson, A., and Schmidt, S. (1971a). Discrete square root filtering: A survey
of current techniques. IEEE Transactions on Automatic Control, AC-16(6): 727–736.
Kaminski, P., Bryson, A., and Schmidt, S. (1971b). Discrete square root filtering: A survey
of current techniques. IEEE Transactions on Automatic Control, AC-16:727–736.
Martin, M. and Kailath, T. (1975). Square-root algorithms for least-squares estimation.
IEEE Transactions on automatic control, AC–20(4): 487–497.
Moon, T. and Stirling, W. (2000). Mathematical Methods and Algorithms for Signal Proces-
sing. Prentice-Hall, Upper Saddle River, New Jersey, USA.
Morf, M., Levy, R., and Kailath, T. (1978). Square-root algorithms for the continuous time
linear least-square estimation problem. IEEE Transactions on Automatic Control, AC-
23(5): 907–911.
Schmidt, S. (1981). The Kalman filter: Its recognition and development for aerospace appli-
cations. Journal of Guidance and Control, 4(1): 4–7.
Verhaegen, M. and Van Dooren, P. (1986). Numerical aspects of different Kalman filter
implementations. IEEE Transactions on Automatic Control, AC-31(l0): 907–917.
Capı́tulo 6
183
184 Capı́tulo 6. Generalización del filtro de Kalman
Puesto que Q es una matriz de covarianza, se sabe que todos sus valores
propios son reales y no negativos. Por tanto, sus valores propios se pueden
denotar como µ2k
Q = DQ̂D⊤ (6.3)
cero, pero en esta sección no se tendrá en cuenta esta suposición. Para tener
más claridad conceptual supóngase, por ejemplo, que el sistema es un avión
en vuelo con el viento golpeándolo. Se utiliza un anemómetro para medir la
velocidad del viento como una entrada al filtro de Kalman. Ası́, las ráfagas
aleatorias de viento afectan tanto al proceso (es decir, la dinámica del avión)
como a la medida (es decir, la velocidad del viento detectada). Se ve que hay
una correlación entre el ruido del proceso y el ruido de las mediciones. De
la ecuación anterior, se observa que el ruido del proceso en el tiempo k está
correlacionado con el ruido de la medición en el tiempo (k + 1); esto es, vk
está correlacionado con wk+1 . Esto es porque vk afecta al estado en el tiempo
(k + 1), ası́ como wk+1 afecta las mediciones en el tiempo (k + 1).
Con el fin encontrar las ecuaciones del filtro de Kalman teniendo en cuenta
el ruido correlacionado del sistema, se definen los errores de estimación de la
siguiente manera:
ϵ− −
k = xk − x̂k
(6.7)
k = xk − x̂k
ϵ+ +
x̂− +
k = Fk−1 x̂k−1 + Gk−1 uk−1
− − (6.8)
k = x̂k + Kk (yk − Hk x̂k )
x̂+
ϵ− −
k = xk − x̂k
= (Fk−1 xk−1 + Gk−1 uk−1 + vk−1 ) − (Fk−1 x̂+
k−1 + Gk−1 uk−1 )
+
= Fk−1 ϵk−1 + vk−1
(6.9)
ϵk = xk − [x̂−
+ −
k + Kk (yk − Hk x̂k )]
= ϵ−
k − Kk (Hk xk + wk − Hk x̂k )
−
= ϵ− −
k − Kk (Hk ϵk + wk )
6.2. Correlación del ruido del proceso y de la medida 187
{ + + ⊤}
k = E ϵk (ϵk )
P+
{[ ][ − ]⊤ }
− − −
= E ϵk − Kk (Hk ϵk + wk ) ϵk − Kk (Hk ϵk + wk )
{( )( − ⊤ )}
= E ϵ− −
k − Kk H k ϵ k − Kk w k (ϵk ) − (ϵ− ⊤ ⊤ ⊤ ⊤ ⊤
k ) Hk Kk − wk Kk
= P− − ⊤ ⊤ − ⊤ ⊤ − − ⊤ ⊤
k − Pk Hk Kk − E {ϵk wk }Kk − Kk Hk Pk + Kk Hk Pk Hk Kk +
+ Kk Hk E {ϵ− ⊤ ⊤ − ⊤ − ⊤ ⊤ ⊤
k wk }Kk − Kk E {wk (ϵk ) } + Kk E {wk (ϵk ) }Hk Kk +
+ Kk E {wk wk⊤ }K⊤
k (6.11)
∂Tr(ABA⊤ )
= 2AB si B es simétrica
∂A
Y además del Apéndice A, las expresiones (A.2) y (A.4) que se reescriben a
continuación
∂ ∂
Tr(XA) = A⊤ ; Tr(AX⊤ ) = A
∂X ∂X
A partir de estos datos se obtiene la siguiente expresión para la derivada
∂Tr(P+
k)
= −2(I − Kk Hk )P− ⊤ ⊤ ⊤
k Hk + 2Kk (Hk Sk + Sk Hk + Rk ) − 2Sk
∂Kk
[ ]
= 2 Kk (Hk P− ⊤ ⊤ ⊤ − ⊤
k Hk + Hk Sk + Sk Hk + Rk ) − Pk Hk − Sk (6.14)
Kk = (P− ⊤ − ⊤ ⊤ ⊤
k Hk + Sk )(Hk Pk Hk + Hk Sk + Sk Hk + Rk )
−1
(6.15)
Esta es la matriz de ganancia de Kalman óptima para el sistema con ruido del
proceso y de la medida correlacionados. La covarianza del error de estimación
se obtiene de la Ecuación (6.13) haciendo Šk = Hk P− ⊤ ⊤ ⊤
k Hk +Hk Sk +Sk Hk +Rk
6.2. Correlación del ruido del proceso y de la medida 189
como
− ⊤ ⊤ ⊤ ⊤
k = (I − Kk Hk )Pk (I − Kk Hk ) + Kk (Hk Sk + Sk Hk + Rk )Kk −
P+
− Sk K⊤ ⊤
k − Kk Sk =
= P− ⊤ − − ⊤ ⊤ ⊤ ⊤
k + Kk Šk Kk − Kk Hk Pk − Pk Hk Kk − Sk Kk − Kk Sk
= P− − ⊤ −1 ⊤ − ⊤ − ⊤ ⊤
k + (Pk Hk + Sk )Šk Šk Kk − Kk (Hk Pk + Sk ) − (Pk Hk + Sk )Kk
= P− − ⊤ ⊤ − ⊤ − ⊤ ⊤
k + (Pk Hk + Sk )Kk − Kk (Hk Pk + Sk ) − (Pk Hk + Sk )Kk
= P− − ⊤
k − Kk (Hk Pk + Sk ) (6.16)
Solución:
Se puede utilizar el método discutido en la Sección 6.1 para simular el ruido
correlacionado. Las ecuaciones del filtro de Kalman dadas anteriormente se pueden
190 Capı́tulo 6. Generalización del filtro de Kalman
0 = E {x0 }
x̂+
{ + ⊤
} (6.18)
P0 = E (x0 − x̂+
+
0 )(x0 − x̂0 )
⊤ − −1 ⊤ − −1 −1
= P+k (Hk + (Pk ) Sk )(Rk − Sk (Pk ) Sk )
x̂− +
k = Fk−1 x̂k−1 + Gk−1 uk−1
− −
k = x̂k + Kk (yk − Hk x̂k )
x̂+
− ⊤ ⊤ ⊤ (6.19)
k = (I − Kk Hk )Pk (I − Kk Hk ) − Sk Kk − Kk Sk +
P+
+ Kk (Hk Sk + S⊤ ⊤ ⊤
k Hk + Rk )Kk
[ − −1
= (Pk ) + (H⊤ − −1 ⊤ − −1
k + (Pk ) Sk )(Rk − Sk (Pk ) Sk )
−1
×
]−1
× (Hk + S⊤ − −1
k (Pk ) )
= P− − ⊤
k − Kk (Hk Pk + Sk )
6.2. Correlación del ruido del proceso y de la medida 191
ejecutar a continuación para obtener una estimación del estado. En la Figura 6.2
se muestra la respuesta a diferentes niveles de correlación, con correlación cero y
con valores de ±0.25.
S = 0, S =0 S = 0, SFiltro = 0
Filtro
0.0915
0.0914
0.0913
0.0912
0.0911
0.091
0.0909
10 20 30 40 50 0 10 20 30 40 50
Tiempo Tiempo
0.024
0.0238
0.0236
0.0234
10 20 30 40 50 0 10 20 30 40 50
Tiempo Tiempo
(b) Ganancia de Kalman y covarianza del error de estimación con correlación (positiva).
S = −0.25, SFiltro = −0.25 S = −0.25, SFiltro = −0.25
0.0655
0.065
0.0645
0.064
0.0635
0.063
0.0625
10 20 30 40 50 0 10 20 30 40 50
Tiempo Tiempo
(c) Ganancia de Kalman y covarianza del error de estimación con correlación (negativa).
Figura 6.1.: Simulación del filtro de Kalman usando ruido del proceso y de
la medida correlacionados. S = correlación entre el ruido del
proceso y el ruido de la medida. SF ilter = valor de S usado en
el filtro de Kalman.
La Tabla 6.1 muestra (para varios valores de S) la varianza del error de estimación
para el filtro de Kalman estándar (cuando se asume S = 0) y para el filtro de Kal-
man con ruido correlacionado (cuando se utiliza el valor correcto de S). Cuando
S = 0, las varianzas de error de estimación son las mismas para los dos filtros,
192 Capı́tulo 6. Generalización del filtro de Kalman
como se esperaba. Sin embargo, cuando S ̸= 0, el filtro que utiliza el valor correcto
de S se desempeña notablemente mejor que el filtro que asume incorrectamente
que S = 0.
donde ζk−1 es ruido blanco de media cero, no correlacionado con vk−1 . En este
caso, se puede ver que la covarianza entre vk y vk−1 es igual a
{ ⊤
} { ⊤ ⊤
}
E vk vk−1 = E ψvk−1 vk−1 + ζk−1 vk−1
= ψQk−1 + 0 (6.22)
El cero surge porque vk−1 es independiente de ζk−1 y ζk−1 tiene media cero.
Se ve que vk es el ruido coloreado del proceso (ya que se correlaciona con sı́
mismo en otros pasos de tiempo). Entonces se pueden combinar las ecuaciones
(6.20) y (6.21) para obtener
[ ] [ ][ ] [ ]
xk F I xk−1 0
= +
vk 0 ψ vk−1 ζk−1
x̌k = F̌x̌k−1 + v̌k−1 (6.23)
Este es un sistema aumentado con un nuevo estado x̌, una nueva matriz del
sistema F̌ y un nuevo vector de ruido del proceso v̌ cuya covarianza está dada
por
[ ]
{ } 0 { 0 }
E v̌k v̌k⊤ = 0 E ζk ζ ⊤
k
= Q̌k (6.24)
Ahora se escribe la señal yk en función de una señal auxiliar y̌k como sigue
De donde,
donde
Se ve que se tiene una nueva ecuación de medición para la medida y̌k−1 , la cual
tiene una matriz de medición Ȟk−1 y un ruido de medida w̌k−1 . Por lo tanto,
el sistema equivalente nuevo puede escribirse como
Kk = P− ⊤ − ⊤
k Ȟk (Ȟk Pk Ȟk + Rk )
−1
Sk = Qk H⊤
k+1 (6.40)
Ck = Sk (Ȟk P− ⊤
k Ȟk + Rk )
−1
− ⊤ ⊤
k = (I − Kk Ȟk )Pk (I − Kk Ȟk ) + Kk Rk Kk
P+
P− + ⊤ ⊤ ⊤ ⊤ ⊤
k+1 = Fk Pk Fk + Qk − Ck Sk − Fk Kk Sk − Sk Kk Fk
x̂−
k = E {xk |y1 , · · · , yk } (6.41)
No se darán los detalles aquı́, pero en [Bryson and Henrikson, 1968], se mues-
tra que esta minimización conduce a las ecuaciones del estimador dadas en el
Algoritmo 6.3.
Ejemplo 6.2
[ ] [ ]
0.70 −0.15 0.15
xk = xk−1 + v
0.03 0.79 0.21 k−1
[ ]
1 0
yk = x + wk
0 1 k
wk = ψwk−1 + ζk−1
6.3. Ruido coloreado de proceso y medición 199
{ }
E vk vj⊤ = 1δk−j
{ } [ ]
⊤ 0.05 0
E ζk ζj = δ
0 0.05 k−j
{ }
E vk ζj⊤ = 0
Tabla 6.2.: Cuando se incrementa el contenido de color del ruido de medición (es
decir, cuando ψ se incrementa), los filtros para ruido coloreado mejoran
su desempeño en comparación con el filtro estándar de Kalman.
Sin embargo, cuando ψ se incrementa (es decir, el color del ruido de medición au-
menta) se ve que los filtros que toman esto en cuenta proporcionan cada vez un
mejor rendimiento en comparación con el filtro estándar de Kalman. En la Fi-
gura 6.2, se puede observar la respuesta en el tiempo de los filtros de Kalman
estándar, aumentado y Bryson y Henrikson, para el caso de ψ = 0.95. Nótese que
hay mayor dispersión en el caso del filtro de Kalman estándar comparado con
los otros casos. Este ejemplo muestra la mejora que es posible en el desempeño
aplicando los filtros para ruido coloreado de la medida.
200 Capı́tulo 6. Generalización del filtro de Kalman
−1
−2
Error de medición Error de medición
Error de estimación Error de estimación
−3
100 200 300 400 500 0 100 200 300 400 500
Tiempo Tiempo
−1
−2
Error de medición Error de medición
Error de estimación Error de estimación
−3
100 200 300 400 500 0 100 200 300 400 500
Tiempo Tiempo
−1
−2
Error de medición Error de medición
Error de estimación Error de estimación
−3
100 200 300 400 500 0 100 200 300 400 500
Tiempo Tiempo
Figura 6.2.: Simulación del filtro de Kalman usando ruido del proceso y de
la medida coloreados.
6.4. Filtrado de estado estacionario 201
x̂− +
k = Fx̂k−1
− −
k = x̂k + K∞ (yk − Hx̂k ) = Fx̂k−1 + K∞ (yk − HFx̂k−1 )
x̂+ + +
= (I − K∞ H)Fx̂+k−1 + K∞ yk (6.46)
P− + ⊤
k+1 = FPk F + Q (6.47)
P− − ⊤ − ⊤ ⊤ ⊤
k+1 = FPk F − FKk HPk F − FKk S F + Q (6.48)
para obtener
P− − ⊤ − ⊤ − ⊤ ⊤ ⊤ −1 − ⊤
k+1 = FPk F − F(Pk H + S)(HPk H + HS + S H + R) HPk F −
− F(P− ⊤ − ⊤ ⊤ ⊤ −1 ⊤ ⊤
k H + S)(HPk H + HS + S H + R) S F + Q
= FP− ⊤ − ⊤ − ⊤ ⊤ ⊤ −1
k F − F(Pk H + S)(HPk H + HS + S H + R) ×
× (HP− ⊤ ⊤
k + S )F + Q (6.49)
Si P− − −
k converge a un valor de estado estacionario, entonces Pk = Pk+1 para k
grande. Se indicará este valor de estado estacionario como P∞ , lo cual significa
que se puede escribir
Hay sistemas para los cuales la ecuación de Riccati (y aquı́ la ganancia de Kal-
man) no converge al valor de estado estacionario. Además, ésta puede convergir
a diferentes valores de estado estacionario dependiendo de la condición inicial
P0 . Finalmente, aun cuando ésta converge a un valor de estado estacionario,
puede resultar en un filtro de Kalman inestable. Estos temas constituyen un
rico campo de estudio que se ha reportado extensamente en muchos libros y
documentos [Anderson and Moore, 1979; Goodwin and Sin, 1984; Chui and
Chen, 1987; Kailath, 1981; McGarty, 1974]. Se resumirán a continuación los
resultados más importantes de convergencia de la ecuación de Riccati, pero
primero se necesita definir qué significa para un sistema ser controlable sobre
el cı́rculo unitario.
Definición 6.1
El par matricial (F, G) es controlable sobre el cı́rculo unitario si existe alguna
matriz K tal que (F − GK) no posee ningún autovalor con magnitud 1.
Ejemplo 6.3
Considérese el sistema escalar
xk+1 = xk (6.52)
Ejemplo 6.4
Considérese el sistema escalar
xk+1 = 2xk (6.53)
Ejemplo 6.5
Considérese el sistema escalar
[ ] [ ]
F1 0 0
xk+1 = x + u (6.54)
0 1 k 1 k
lazo cerrado tengan magnitud no unitaria y el sistema es, por lo tanto, controlable
sobre el circulo unitario.
k = (I − K∞ H)Fx̂k−1 + K∞ yk
x̂+ +
(6.58)
Más aún, exactamente una de las soluciones DARE definidas positivas resultan
en un filtro de Kalman de estado estacionario estable.
Teorema 6.4
La DARE tiene al menos una solución semidefinida positiva P∞ si y solo si (F, H)
es detectable. Además, al menos una de tales soluciones resulta en un filtro de
Kalman de estado estacionario marginalmente estable.
P = F P F ⊤ − F P H ⊤ (HP H ⊤ + R)−1 HP F ⊤ + Q
= P − P (P + 1)−1 P + 1
De aquı́ se obtiene
(P + 1)−1 P = 1, es decir, P 2 + P − 1 = 0
Resolviendo la ecuación cuadrática anterior se obtiene
√
1± 5
P =
2
Ası́, la DARE tiene dos soluciones, una negativa y la otra positiva. Si se utiliza la
solución negativa en el filtro de Kalman de estado estacionario se obtiene
√
⊤ ⊤ −1 1− 5
K = P H (HP H + R) = √
3− 5
2
k = (I − KH)F x̂k−1 + Kyk =
x̂+ +
√ x̂+k−1 + Kyk ≈ 2.62x̂k−1 + Kyk
+
3− 5
208 Capı́tulo 6. Generalización del filtro de Kalman
K = P H ⊤ (HP H ⊤ + R)−1
√
1+ 5
= √
3+ 5
2
k = (I − KH)F x̂k−1 + Kyk =
x̂+ +
√ x̂+ + Kyk
3 + 5 k−1
≈ 0.382x̂+
k−1 + 0.618yk
Ejemplo 6.7
Considérese un sistema escalar con F = 1, H = 1, Q = 0, R = 1 y S = 0.
Nótese que (F, H) es detectable. Sin embargo, (F, G) no es controlable en el cı́rculo
unitario para todo G tal que GG⊤ = Q. Se sabe del Teorema 6.2 que la DARE
no tiene una solución semidefinida positiva que resulte en un filtro de Kalman
estable. Sin embargo, se sabe del Teorema 6.4 que la DARE tiene una solución
semidefinida positiva que resulta en un filtro de Kalman marginalmente estable.
La DARE para este sistema está dada por
P = F P F ⊤ − F P H ⊤ (HP H ⊤ + R)−1 HP F ⊤ + Q
= P − P (P + 1)−1 P
Ésta tiene dos soluciones para P , ambas 0 (es decir, semidefinida positiva). Si se
usa esta solución en el filtro de Kalman en estado estacionario se obtiene
K=0
x̂+ +
k = x̂k−1
Ejemplo 6.8
Considérese un sistema escalar con F = 2, H = 1, Q = 0, R = 1 y S = 0. Nótese
que (F, H) es detectable. También (F, G) es controlable sobre y en el interior del
circulo unitario para todo G tal que GG⊤ = Q. Se sabe del Teorema 6.2 que la
DARE tiene exactamente una solución semidefinida positiva que resulta en un
filtro de Kalman estable.
Sin embargo, se sabe del Teorema 6.4 que la DARE tiene una solución semi-
definida positiva que resulta en un filtro de Kalman marginalmente estable que
es estable. También se sabe del Teorema 6.3 que esta solución DARE es definida
6.4. Filtrado de estado estacionario 209
P = F P F ⊤ − F P H ⊤ (HP H ⊤ + R)−1 HP F ⊤ + Q
= 4P − 4P (P + 1)−1 P
es decir,
4P 2
3P − =0
P +1
P 2 − 3P = 0
o sea,
P =0 y P =3
Como se puede observar hay dos soluciones para P , una de las cuales es 0 (es
decir, semidefinida positiva) y la otra es 3 (es decir, definida positiva). Si se toma
P = 0, en el filtro de Kalman de estado estacionario, se obtiene
K=0
x̂+ +
k = 2x̂k−1
Factorizando F y F⊤ del principio y del final de los tres primeros términos del
lado derecho conduce a
[
Pk+1 = F Pk − Pk H⊤ R−1 HPk +
] ⊤
+ Pk H⊤ R−1 H(H⊤ R−1 H + P−1 )−1 ⊤ −1
H R HP k F +Q
{ [ k
]}
= F Pk − Pk H R H Pk − (H R H + Pk ) H R HPk F⊤ + Q
⊤ −1 ⊤ −1 −1 −1 ⊤ −1
{
= F Pk − Pk H⊤ R−1 H(H⊤ R−1 H + P−1 k )
−1
×
[ ⊤ −1 ]} ⊤
× (H R H + P−1 ⊤ −1
k )Pk − H R HPk F +Q
{ ⊤ −1 ⊤ −1 −1 −1
} ⊤
= F Pk − Pk H R H(H R H + Pk ) F +Q
[ ⊤ −1 ] ⊤ −1
= FPk (H R H + Pk ) − H R H (H R H + P−1
−1 ⊤ −1 −1 ⊤
k ) F +Q
= F(H⊤ R−1 H + P−1 −1 ⊤
k ) F +Q
= FPk (H⊤ R−1 HPk + I)−1 F⊤ + QF−⊤ F⊤
[ ]
= FPk + QF−⊤ (H⊤ R−1 HPk + I) (H⊤ R−1 HPk + I)−1 F⊤
[ ]
= (F + QF−⊤ H⊤ R−1 H)Pk + QF−⊤ (H⊤ R−1 HPk + I)−1 F⊤ (6.62)
Pk = Υk Z−1
k (6.63)
6.4. Filtrado de estado estacionario 211
Estas ecuaciones para Υk+1 y Zk+1 se pueden escribir como una ecuación única:
[ ] [ −⊤ ][ ]
Zk+1 F F−⊤ H⊤ R−1 H Zk
=
Υk+1 QF−⊤ F + QF−⊤ H⊤ R−1 H Υk
[ ]
Z
=H k (6.66)
Υk
Yk+1 = DYk
[ ] [ −1 ][ ]
Y1,k+1 Λ 0 Y1k
= (6.72)
Y2,k+1 0 Λ Y2k
Y2,k+1 = ΛY2k
(6.73)
Y2k = Λk Y2,0
Similarmente se ve que
Υk = Ψ22 Ψ−1
12 Zk (6.80)
214 Capı́tulo 6. Generalización del filtro de Kalman
P− −1
∞ = Ψ22 Ψ12 (6.79)
◃ Nótese que Ψ12 debe ser invertible para que este método funcione.
Υ k = P k Zk (6.81)
Ejemplo 6.10
Considérese el sistema escalar de la Ecuación (6.44):
xk+1 = xk + vk vk ∼ N (0, 1)
yk = xk + wk wk ∼ N (0, 1)
Se ve que F = H = Q = R = 1. Sustituyendo esos valores en la expresión para la
matriz Hamiltoniana da
[ −⊤ ] [ ]
F F −⊤ H ⊤ R−1 H 1 1
H= =
QF −⊤ F + QF −⊤ H ⊤ R−1 H 1 2
Los valores propios de H son 0.38 y 2.62. Ninguno de los valores propios tiene
magnitud uno por lo que se está en condiciones de continuar con el procedimiento.
El autovector de H que corresponde al valor propio fuera del cı́rculo unitario es
⊤
[0.5257 0.8507] . Se forma la correspondiente matriz de autovectores como
[ ] [ ]
Ψ12 0.5257
=
Ψ22 0.8507
Nótese que Ψ12 es invertible por lo que se podrá continuar con el problema. La
solución de la ecuación de Riccati de estado estacionario es
6.5. Resumen 217
0.8507
P = Ψ22 Ψ−1
12 = = 1.62
0.5257
La ganancia de Kalman de estado estacionario por lo tanto se calcula de la Ecua-
ción (6.19) como
(1.62)(1)
K = P H ⊤ (HP H ⊤ + R)−1 = = 0.62
(1)(1.62)(1) + 1
la cual coincide con la Ecuación (6.45).
6.5. Resumen
En este capı́tulo, se discutió una variedad de generalizaciones del filtro de Kal-
man que hacen que el filtro sea aplicable ampliamente a un número mayor
de problemas. Ruidos del proceso y de la medida correlacionados y coloreados
fueron estudiados desde los primeros dı́as del invento del filtro de Kalman. En
este capı́tulo se ha demostrado que las modificaciones del filtro teniendo en
cuenta la correlación y color puede mejorar el rendimiento de la estimación.
Sin embargo, si estos métodos son dignos de la complejidad adicional y esfuer-
zo computacional depende del problema. Una de las extensiones más prácticas
del filtro Kalman lo constituye el filtro de Kalman de estado estacionario. El
filtro de Kalman de estado estacionario a menudo se desempeña casi de la mis-
ma forma como el filtro variable en el tiempo teóricamente más riguroso. Sin
embargo, el filtro de estado estacionario requiere sólo una fracción del costo
computacional. Para el análisis del filtro de estado estacionario se aplica el
modelo Hamiltoniano el cual permite determinar la solución en estado esta-
cionario de la ecuación de Riccati. Este método tiene algunas limitaciones, es
decir, si la matriz Hamiltoniana tiene valores propios con magnitud igual a
uno, entonces el análisis falla.
Bibliografı́a
Anderson, B. and Moore, J. (1979). Optimal Control. Information and System Sciences
Series. Prentice–Hall, Inc., Englewood Cliffs, New Jersey 07632, USA.
Arnold III, W. F. and Laub, A. J. (1984). Generalized eigenproblem algorithms and software
for algebraic Riccati equations. Proceedings of the IEEE, 72(12): 1746–1754.
Bitmead, R., Gevers, M., Petersen, I., and Kaye, R. (1985). Monotonicity and stabiliza-
bility properties of solutions of the Riccati difference equation: Propositions, lemmas,
theorems, fallacious conjectures and counterexamples. Systems and Control Letters,
5(5): 309–315.
Bryson, A. and Henrikson, L. (1968). Estimation using sampled data containing sequentially
correlated noise. Journal of Spacecraft and Rockets, 5(6): 662–665.
Bryson, A. and Johanson, D. (1965). Linear filtering for time-varying systems using mea-
surements containing colored noise. IEEE Transactions on Automatic Control, AC-
10(1): 4–10.
Bucy, R. and Joseph, P. (1968). Filtering for Stochastic Processes with Applications to
Guidance. Number 23 in Interscience tracts in pure and applied mathematics. John
Wiley & Sons, New York, NY, USA.
Chui, C. and Chen, G. (1987). Kalman Filtering with Real-Time Applications. Springer-
Verlag, New York, USA.
Crassidis, J. L. and Junkins, J. L. (2012). Optimal Estimation of Dynamic Systems. 6000
Broken Sound Parkway NW, Suite 300 Boca Raton, FL 33487-2742, USA., 2a edition.
Goodwin, G. and Sin, K. (1984). Adaptive Filtering, Prediction and Control. Prentice Hall,
Englewood Cliffs, New Jersey, USA.
Kailath, T. (1981). Lectures on Wiener and Kalman Filtering. Springer–Verlag, New York,
USA.
Kailath, T., Sayed, A., and Hassibi, B. (2000). Linear Estimation. Prentice-Hall, Upper
Saddle River, New Jersey, USA.
Kalman, R. E. (1960). A new approach to linear filtering and prediction problems. Tran-
sactions of the ASME–Journal of Basic Engineering, 82(Series D): 35–45.
Maybeck, P. (1979). Stochastic Models, Estimation, and Control, volume Volume 1. Acade-
mic Press, New York.
219
220 BIBLIOGRAFı́A
McGarty, T. (1974). Stochastic Systems and State Estimation. John Wiley & Sons, New
York, USA.
Poubelle, M., Petersen, I., Gevers, M., and Bitmead, R. (1986). A miscellany of results
on an equation of Count J. F. Riccati. IEEE Transactions on Automatic Control,
AC-31(7): 651–654.
Simon, D. (2006). Optimal State Estimation: Kalman, H-infinity, and Nonlinear Approaches.
John Wiley & Sons, Inc., New Jersey, USA.
Stear, E. and Stubberud, A. (1968). Optimal filtering for Gauss–Markov noise. International
Journal of Control, 8(2): 123–130.
Stengel, R. (1994a). Optimal Control and Estimation. Dover Publications, Inc., New York.
Stengel, R. (1994b). Optimal Control and Estimation. New York, USA.
Vaughan, D. (1970). A nonrecursive algebraic solution for the discrete Riccati equation.
IEEE Transactions on Automatic Control, 15(5): 597 – 599.
Capı́tulo 7
221
222 Capı́tulo 7. El filtro de Kalman continuo
xk = xk−1 + vk−1
vk ∼ (0, Q) (7.1)
x0 = 0
xk = v0 + v1 + · · · + vk−1 (7.2)
Supóngase que se tiene una medida en tiempo discreto de una constante x cada
T segundos. Los tiempos de medida son tk = kT (k = 1, 2, . . .):
xk = xk−1
yk = xk + wk (7.12)
wk ∼ (0, R)
+ Pk+ R
Pk+1 = (7.13)
Pk+ + R
7.1. Ruido blanco en tiempo discreto y en tiempo continuo 225
wk ∼ (0, R)
(7.17)
w(t) ∼ (0, Rc )
Los resultados de las secciones anteriores se pueden combinar con los resultados
de la Sección 2.4 para obtener una simulación discretizada de un sistema de
tiempo continuo ruidoso, con el propósito de implementar un estimador de
226 Capı́tulo 7. El filtro de Kalman continuo
ẋ = Ax + Bu + v
y = Cx + w
(7.19)
v ∼ (0, Qc )
w ∼ (0, Rc )
Tanto v(t) y w(t) son ruidos en tiempo continuo y u(t) es una entrada co-
nocida. Este sistema es equivalente aproximadamente al siguiente sistema en
tiempo discreto
∫ T
AT
xk = e xk−1 + e AT
e−Aτ dτ Buk−1 + vk−1
[0 ]
AT
= e xk−1 + e AT
I − e−AT A−1 Buk−1 + vk−1
(7.20)
yk = Cxk + wk
vk ∼ (0, Qc T)
wk ∼ (0, Rc /T)
La ganancia del filtro de Kalman de tiempo discreto para este sistema se dedujo
en la Sección 4.4 como
Kk = P−1 ⊤ − ⊤
k H (HPk H + R)
−1
(7.24)
P− + ⊤
k+1 = (I + AT)Pk (I + AT) + Qc T
+ ⊤ + ⊤ 2
= P+ +
k + (APk + Pk A + Qc )T + APk A T
Sustrayendo P−
k de ambos lados y dividiendo entre T da
P− −
k+1 − Pk −Kk CP− ⊤ − − ⊤ AKk CP−
= k
+ AP+
k A T + APk + P k A + Qc − k
T
T T T
Kk CP−kA
⊤
− T
T
donde no se han simplificado los dos últimos términos de la derecha para efectos
de tomar el lı́mite.
Tomando el lı́mite cuando T → 0 y usando la Ecuación (7.25) da
P− −
k+1 − Pk
Ṗ = lı́m = −PC⊤ R−1 ⊤
c CP + AP + PA + Qc
(7.27)
T→0 T
Esta expresión se denomina ecuación diferencial de Riccati en P y puede usarse
para calcular la covarianza del error de estimación del filtro de Kalman de
tiempo continuo. Esto requiere n2 integraciones porque P es una matriz n × n.
Pero P es simétrica, ası́ en la práctica sólo se necesita integrar n(n + 1)/2
ecuaciones con el fin de resolver P. En la Sección 4.4 se dedujeron las ecuaciones
del filtro de Kalman para x̂ como
x̂− +
k = Fx̂k−1 + Guk−1
− − (7.28)
k = x̂k + Kk (yk − Hx̂k )
x̂+
x̂+ + +
k = x̂k−1 + ATx̂k−1 + BTuk−1 +
[ ]
PC⊤ R−1c T yk − Cx̂k−1 − CATx̂k−1 − CBTuk−1
+ +
(7.30)
Substrayendo x̂+
k−1 en ambos lados, dividiendo por T, y tomando el lı́mite
cuando T → 0, da
T→0 T
= Ax̂ + Bu + PC⊤ R−1
c (y − Cx̂) (7.31)
Esta expresión corresponde a la ecuación diferencial que puede ser usada para
integrar el estado estimado en el filtro de Kalman de tiempo continuo. El
filtro de Kalman en tiempo continuo puede resumirse como se muestra en el
Algoritmo 7.1.
Ejemplo 7.1
En este ejemplo se usará el filtro de Kalman de tiempo continuo para estimar una
constante, dadas las medidas de ruido en tiempo continuo:
ẋ = 0
y =x+w
w ∼ (0, R)
ẋ = Ax + Bu + v
y = Cx + w
(7.33)
v ∼ (0, Qc )
w ∼ (0, Rc )
x̂(0) = E {x(0)}
{ }
P(0) = E (x(0) − x̂(0))(x(0) − x̂(0))⊤
K = PC⊤ R−1 c (7.34)
˙x̂ = Ax̂ + Bu + K(y − Cx̂)
Ṗ = −PC⊤ R−1 ⊤
c CP + AP + PA + Qc
Entonces
Ṗ = −P 2 /R
con la condición inicial P (0) = P0 . De esto se puede obtener
dP −dτ
=
P2 R
∫ P (t) ∫ t
1 1
2
dP = − dτ
P (0) P 0 R
( ) t
− P −1 − P0−1 = −
( R )
−1 −1 t
P = P0 +
R
( )
−1 t −1
P = P0 +
R
P0
=
1 + P0 t/R
7.1. Ruido blanco en tiempo discreto y en tiempo continuo 231
lı́m P = 0
t→∞
K = P C ⊤ R−1
P0 /R
=
1 + P0 t/R
lı́m K = 0
t→∞
Esto muestra que cuando el tiempo tiende a infinito, x̂ alcanza un valor de esta-
do estacionario. Esto es intuitivo porque cuando se tiene un número infinito de
mediciones de una constante, el estimado de esa constante llega a ser perfecto y
medidas adicionales no pueden mejorar el estimado. Además, la ganancia de Kal-
man tiende a cero cuando el tiempo tiende a infinito, lo que a su vez dice que se
tienen que ignorar mediciones adicionales (ya que el estimado llega a ser perfecto).
Finalmente, la covarianza P tiende a cero cuando el tiempo tiende a infinito, lo
cual dice que la incertidumbre en el estimado tiende a cero, significando que el
estimado es perfecto.
Ejemplo 7.2
En este ejemplo se tienen las medidas de la velocidad de un objeto que se mueve en
una dimensión. El objeto está sujeto a aceleraciones aleatorias. Se desea estimar la
velocidad x a partir de medidas ruidosas de la misma. Las ecuaciones del sistema
y de la medida están dadas como
ẋ = v
y =x+w
v ∼ (0, Q)
w ∼ (0, R)
232 Capı́tulo 7. El filtro de Kalman continuo
Solución:
Ṗ = −P C ⊤ R−1 CP + AP + P A⊤ + Q
= −P 2 /R + Q
con la condición inicial P (0) = P0 . De esto se puede obtener
( )
dP = Q − P 2 /R dτ
∫ t ∫ P (t)
dP
dτ = (7.36)
P (0) Q − P /R
2
0
∫ t ∫ P (t)
dP
dτ = √ √ √ √
0 P (0) ( Q + P/ R)( Q − P/ R)
∫ t [∫ ∫ P (t) ]
P (t)
1 dP dP
dτ = √ √ √ + √ √
0 2 Q P (0) Q + P/ R P (0) Q − P/ R
√ (√ √ )P (t)
R Q + P/ R
t = √ ln √ √
2 Q Q − P/ R P (0)
√ √ √
R ( QR + P )( QR − P0 )
t = √ ln √ √
2 Q ( QR − P )( QR + P0 )
√ √ √
( QR + P )( QR − P0 )
e 2 Q/Rt
= √ √
( QR − P )( QR + P0 )
Resolviendo esto para P da
[ √ √ √ ]
√ P0 − QR + ( QR + P0 )e2 Q/Rt
P = QR √ √ √
QR − P0 + ( QR + P0 )e2 Q/Rt
√
lı́m P = QR
t→∞
Por lo tanto,
( ) ( ) √
−1 P −1 P Q
tanh √ = tanh √ 0 + t
QR QR R
Resolviendo para P
[ ( ) √ ]
√ −1 P Q
P = QR tanh tanh √ 0 + t
QR R
Λ(0) = P(0)
(7.45)
Y(0) = I
donde las matrices ϕij están definidas como las cuatro submatrices n × n en
eJT . De la factorización original se supone que se tiene Λ = PY, por lo tanto
esta ecuación se puede escribir como
[ ] [ ][ ]
Λ(t + T) ϕ11 (T) ϕ12 (T) P(t)Y(t)
= (7.48)
Y(t + T) ϕ21 (T) ϕ22 (T) Y(t)
Puesto Λ(t + T) = P(t + T)Y(t + T), la primera ecuación se puede escribir como
P(t + T) = [ϕ11 (T)P(t) + ϕ12 (T)] [ϕ21 (T)P(t) + ϕ22 (T)]−1 (7.52)
Esta puede ser una forma rápida para resolver P en lugar de integrar la ecua-
ción de Riccati. Se debe tener en cuenta que con este método no hay que
preocuparse por el tamaño del paso de integración. Este método puede usarse
para propagar desde P(t) hasta P(t + T) en una sola ecuación, para cualquier
valor de t y de T.
Ejemplo 7.3
Supóngase que se desea estimar la tasa desplazamiento ξ de un giroscopio (se
asume que es constante) dadas las medidas del ángulo de giro θ. El modelo del
7.2. Soluciones alternativas para la ecuación de Riccati 237
θ̇ = ξ
y =θ+w
[ ] [ ][ ]
θ̇ 0 1 θ
= (7.53)
ξ˙ 0 0 ξ
[ ]
θ
y = [1 0] +w
ξ
w ∼ (0, R)
Solución:
[ ] [ ]
Λ̇ Λ
=J
Ẏ Y
[ ]
A Q
J=
C⊤ R−1 C −A⊤
0 1 0 0
0 0 0 0
=
1/R 0 0 0
0 0 −1 0
La matriz de transición para la ecuación diferencial de Λ y Y se calcula como
1 t 0 0
0 1 0 0
e(Jt) =
t/R t2 /2R 1 0
−t2 /2R −t3 /6R −t 1
[ ]
ϕ11 (t) ϕ12 (t)
=
ϕ21 (t) ϕ22 (t)
Esto es, la incertidumbre tiende a cero cuando el tiempo tiende a infinito. Es-
to ocurre por que el ruido del proceso es cero (es decir, se está estimando una
constante). Puesto que K = PC⊤ R−1 , se ve que la ganancia de Kalman también
tiende a cero cuando el tiempo tiende a infinito. Esto simplemente significa que
eventualmente se han tomado tantas medidas que el conocimiento es completo.
Medidas adicionales no pueden dar ninguna información nueva, por lo tanto se
ignoran las nuevas medidas.
Γ−1 Γ̇ = MS
(7.56)
Γ̇⊤ Γ−⊤ = MI
Γ̇ = ΓMS (7.57)
Esta es llamada una ecuación algebraica de Riccati (ARE). Para ser más es-
pecı́fico, esta se denomina ARE continua (CARE).
La solución de la ARE puede no siempre existir, e incluso si existe puede no
resultar en un filtro de Kalman estable. Se resumirá el resultado más impor-
tante de convergencia de la ecuación de Riccati a continuación, pero primero
se necesita definir qué significa para un sistema ser controlable sobre el eje
imaginario.
Definición 7.1
El par de matrices (A, B) es controlable sobre el eje imaginario si existe alguna
matriz K tal que (A − BK) no tiene autovalores sobre el eje imaginario.
Se supone que Q ≥ 0 y R > 0. Se define G como una matriz cualquiera tal que
GG⊤ = Q. La ganancia de Kalman de estado estacionario correspondiente K
se da como
Este teorema es análogo al Teorema 6.1 para filtros de Kalman de tiempo dis-
creto. El teorema anterior no excluye la existencia de soluciones CARE que son
definidas negativas o indefinida. Si existen tales soluciones, entonces resultarı́a
en un filtro de Kalman inestable. Si se debilita la condición de estabilizabilidad
en el teorema anterior, se obtiene el siguiente
Teorema 7.2
La CARE tiene al menos una solución P semidefinida positiva si y solo si las
siguientes dos condiciones se cumplen.
1. (A, C) es detectable
2. (A, G) es controlable sobre el eje imaginario.
Más aún, exactamente una de las soluciones ARE semidefinida positiva resulta en
un filtro de Kalman de estado estacionario estable.
Este teorema es análogo al Teorema 6.3 para el filtro de Kalman de tiempo dis-
creto. Si se deja la condición de controlabilidad en los dos teoremas anteriores,
se obtiene el siguiente.
Teorema 7.4
La CARE tiene al menos una solución semidefinida positiva P si (A, C) es de-
tectable. Además, al menos una de tales soluciones resulta en un filtro de Kalman
de estado estacionario marginalmente estable.
Este teorema es análogo al Teorema 6.4 para filtros de Kalman de tiempo dis-
creto. Nótese que el filtro resultante es solamente estable marginalmente, por lo
tanto, puede tener autovalores sobre el eje imaginario. También nótese que este
teorema posee una condición suficiente (pero no necesaria). Esto es, puede ha-
ber un filtro de Kalman de estado estacionario estable aún sı́ las condiciones del
teorema anterior no se mantienen. Además, aún sı́ las condiciones del teorema
se mantienen, puede haber soluciones CARE que resultan en filtros de Kalman
inestables. Resultados adicionales relacionados al filtro de tiempo continuo de
estado estacionario se pueden encontrar en muchos lugares, incluyendo [Aoki,
1967], [Bucy and Joseph, 1968], [Kwakernaak and Sivan, 1972]. Muchos de los
filtros de Kalman prácticos se aplican a sistemas que no reúnen las condiciones
de los teoremas anteriores, pero los filtros aún en esa condiciones funcionan
bien en la práctica.
Ejemplo 7.4
Se considera el siguiente sistema de dos estados ([Bucy and Joseph, 1968]):
[ ]
a1 0
ẋ = x+v (7.64)
0 a2
7.3. Filtro de Kalman en tiempo continuo de estado estacionario 245
[ ]
1 0
y= x+w
0 1
[ ]
q11 q12
Q= (7.65)
q12 q22
[ ]
r1 0
R=
0 r2
En el resto de este ejemplo, se usa el sı́mbolo G para nombrar cualquier matriz
tal que GG⊤ = Q.
La ecuación diferencial de Riccati para el filtro de Kalman está dada como
Ésta puede ser escrita como las siguientes tres ecuaciones diferenciales acopladas.
q11
γ1 = + a21
r1
q22
γ2 = + a22
r2
esto resulta en un filtro de Kalman de estado estacionario estable.
Si a1 = a2 < 0 y q12 ̸= 0, entonces (A, C) es detectable y (A, G) es estabilizable.
Los resultados del Teorema 7.1, por lo tanto, también se aplican a esta situación. Se
puede demostrar [Bucy and Joseph, 1968], que la única solución ARE semidefinida
246 Capı́tulo 7. El filtro de Kalman continuo
ẋ = −x + v
(7.73)
y = x+w
donde v y w son ruidos blancos del proceso y de la medida de media cero y no
correlacionados, con varianzas respectivas de Q = 2 y R = 1. El filtro de Kalman
de estado estacionario para este sistema se puede obtener resolviendo la Ecuación
(7.34) con Ṗ = 0, de la cual se obtiene
√ √
x̂˙ = − 3x̂ + ( 3 − 1)y (7.74)
Tomando la transformada de Laplace de este estimador da
√ √
(s + 3)X̂(s) = ( 3 − 1)Y (s) (7.75)
En otras palabras, el filtro de Kalman es equivalente a pasar la medición y(t) a
través de la función de transferencia G(s), lo cual está dado como
√
3−1
G(s) = √ (7.76)
s+ 3
La respuesta al impulso del filtro de Kalman es obtenida tomando la inversa de
la transformada de la Laplace, lo cual da
√ √
g(t) = ( 3 − 1)e− 3t , t≥0 (7.77)
Ahora se va a obtener el espectro de potencia del estado tomando la trasformada
248 Capı́tulo 7. El filtro de Kalman continuo
7.4. Dualidad
Es interesante observar la dualidad entre la estimación óptima y un control
óptimo. El problema de estimación óptima comienza con las ecuaciones del
sistema y la medición
ẋ = Ax + v
v ∼ N(0, Q)
(7.80)
y = Cx + w
w ∼ N(0, R)
costo
∫ tf { }
Je = E (x − x̂)⊤ (x − x̂) dt (7.81)
0
ẋ = Ax + Cu (7.83)
Pc (tf ) = ϕ(tf )
Ṗc = −A⊤ Pc − Pc A + Pc CR−1 C⊤ Pc − Q
(7.85)
Kc = R−1 C⊤ Pc
u = Kc x
de Riccati tienen la misma forma, excepto que son negativos entre sı́ y A y C
se sustituyen por sus transpuestas. La ganancia estimador Ke , y la ganancia
del controlador Kc tienen formas muy similares. Las matrices de covarianza Q
y R en el problema de estimación tienen duales en las matrices de ponderación
de la función de costo del problema de control óptimo.
La relación dual entre los problemas de estimación y de control se observó en
los primeros trabajos sobre el filtro de Kalman [Kalman, 1960], [Kalman and
Bucy, 1961]. Desde entonces, se ha utilizado muchas veces para extrapolar los
resultados conocidos de un problema con el fin de obtener nuevos resultados
para el problema dual.
7.5. Resumen
En este capı́tulo, se ha deducido el filtro de Kalman de tiempo continuo me-
diante la aplicación de un argumento limitador del filtro de Kalman de tiempo
discreto. Sin embargo, ası́ como hay varias maneras para obtener el filtro de
Kalman tiempo discreto, también hay varias maneras para obtener el filtro de
Kalman de tiempo continuo. La derivación original realizada por Kalman y
Bucy [Kalman and Bucy, 1961] implicaba la solución de la ecuación integral
de Wiener-Hopf.
Se ha visto que tanto la ecuación diferencial como la ecuación algebraica de
Riccati son la clave para la solución del filtro de Kalman de tiempo continuo.
La versión escalar de lo que hoy se conoce como la ecuación de Riccati se
estudió inicialmente por matemáticos como James Bernoulli y John Bernoulli
en los años 1600 y Jacopo Riccati, Daniel Bernoulli, Leonard Euler, Jean-
Lerond d’Alembert y Adrien Legendre en el siglo XVIII. La ecuación fue por
primera vez llamada “la ecuación de Riccati” por d’Alembert en 1763 [Watson,
1944].
El filtro de Kalman en tiempo continuo se aplica a los sistemas con ruido
blanco de tiempo continuo tanto en la ecuación del proceso como en la de
medición. El ruido blanco de tiempo continuo es no intuitivo ya que tiene una
correlación infinita consigo mismo en el instante actual, pero una correlación
cero consigo mismo cuando se separan por tiempos arbitrariamente pequeños
pero distintos de cero. Sin embargo, el ruido blanco de tiempo continuo es un
caso lı́mite de ruido blanco de tiempo discreto, que es intuitivamente acepta-
ble. Por lo tanto, el ruido blanco de tiempo continuo puede ser aceptado como
una aproximación a la realidad. Esto corresponde a muchas otras aproxima-
7.5. Resumen 251
253
254 BIBLIOGRAFı́A
Schmidt, S. (1981). The Kalman filter: Its recognition and development for aerospace appli-
cations. Journal of Guidance and Control, 4(1): 4–7.
Shanahan, W. (1982). Circuit models for prediction, Wiener filtering, Levinson and Kalman
filters for ARMA time series. In IEEE International Conference on Acoustics, Speech,
and Signal Processing ICASSP’82, volume 7, pages 266 – 269.
Simon, D. (2006). Optimal State Estimation: Kalman, H-infinity, and Nonlinear Approaches.
John Wiley & Sons, Inc., New Jersey, USA.
Stengel, R. (1994). Optimal Control and Estimation. New York, USA.
Strang, G. (2006). Linear Algebra and its Applications. Thomson Corporation, Belmont,
CA, USA, 4a edition.
Watson, G. (1944). A Treatise on the Theory of Bessel functions. Cambridge University
Press, Cambridge, United Kingdom, 2a edition.
Capı́tulo 8
255
256 Capı́tulo 8. El filtro de Kalman no lineal
ẋ = f (x, u, v, t)
y = h(x, w, t)
(8.1)
v ∼ (0, Q)
w ∼ (0, R)
∂h ∂h
y ≈ h(x0 , w0 , t) + (x − x0 ) + (w − w0 )
∂x 0 ∂w 0
= h(x0 , w0 , t) + C∆x + M∆w (8.3)
ẋ0 = f (x0 , u0 , v0 , t)
(8.4)
y0 = h(x0 , w0 , t)
∆ẋ = ẋ − ẋ0
(8.5)
∆y = y − y0
258 Capı́tulo 8. El filtro de Kalman no lineal
ẋ0 = f (x0 , u0 , 0, t)
(8.7)
y0 = h(x0 , 0, t)
0 = E {x0 }
x̂+
{ } (8.8)
+ ⊤
P+0 = E (x 0 − x̂+
0 )(x0 − x̂ 0 )
∆y = y − y0 (8.11)
∆x̂(0) = 0
{ }
P(0) = E [∆x(0) − ∆x̂(0)] [∆x(0) − ∆x̂(0)]⊤
∆x̂˙ = A∆x̂ + K(∆y − C∆x̂) (8.12)
K = PC⊤ R̃−1
Ṗ = AP + PA⊤ + Q̃ − PC⊤ R̃−1 CP
x̂ = x0 + ∆x̂ (8.13)
8.2. Filtro de Kalman extendido 259
∆ẋ = A∆x + Lv
= A∆x + ṽ (8.14)
ṽ ∼ (0, Q̃), Q̃ = LQL⊤
∆y = C∆x + Mw
= C∆x + w̃ (8.15)
⊤
w̃ ∼ (0, R̃), R̃ = MRM
Las ecuaciones anteriores corresponden a un sistema lineal con estados ∆x y
medidas ∆y, de manera que se puede utilizar el filtro de Kalman para estimar
∆x. La entrada al filtro consiste en ∆y, la cual es la diferencia entre el estado
real y y el estado nominal y0 . La salida del filtro de Kalman es ∆x, la cual es
un estimado de la diferencia entre el estado real x y el estado nominal x0 . Las
ecuaciones modificadas para el filtro de Kalman linealizado son
∆x̂(0) = 0
{ }
P(0) = E [∆x(0) − ∆x̂(0)] [∆x(0) − ∆x̂(0)]⊤
∆x̂˙ = A∆x̂ + K(∆y − C∆x̂) (8.16)
⊤ −1
K = PC R̃
Ṗ = AP + PA⊤ + Q̃ − PC⊤ R̃−1 CP
x̂ = x0 + ∆x̂
−R ωλ ua + q1
i̇a = ia + sen θ +
L L L
−R ωλ ub + q2
i̇b = ib − cos θ +
L L L (8.21)
−3λ 3λ fω
ω̇ = ia sen θ + ib cos θ − + q3
2J 2J J
θ̇ = ω
donde ia e ib son las corrientes en los dos devanados, θ y ω son la posición angular
y la velocidad del rotor, R y L son la resistencia y la inductancia del devanado,
λ es la constante de flujo y f es el coeficiente de fricción viscosa. Las entradas de
control ua y ub consisten de las tensiones aplicadas a través de los dos devanados
y J es el momento de inercia del eje del motor y la carga. El estado se define como
⊤ ⊤
x = [x1 x2 x3 x4 ] = [ia ib ω θ] (8.22)
Los términos qi son el ruido del proceso debido a incertidumbre en las entradas
de control (q1 y q2 ) y el torque de la carga (q3 ).
Los datos de placa del motor están dados por: resistencia del devanado, R =
1.9 Ω, inductancia del devanado L = 0.003 H, constante del motor, λ = 0.1,
momento de inercia, J = 0.00018 N·m, coeficiente de fricción, f = 0.001 N·s·m−2 ,
I0 = 0.3708 A. Se tomó un paso de integración de dt = 0.0005 s.
8.2. Filtro de Kalman extendido 261
Q̃ = LQL⊤
(8.19)
R̃ = MRM⊤
x̂(0) = E {x(0)}
{ }
P(0) = E [x(0) − x̂(0)] [x(0) − x̂(0)]⊤
x̂˙ = f (x̂, u, v0 , t) + K[y − h(x̂, w0 , t)] (8.20)
K = PC⊤ R̃−1
Ṗ = AP + PA⊤ + Q̃ − PC⊤ R̃−1 CP
donde w(1) y w(2) son procesos independientes de ruido blanco con media cero
y desviaciones estándar asumidas de 100 mA. Las entradas de control nominal se
establecen en
ua (t) = sen(2πt)
(8.25)
ub (t) = cos(2πt)
0.4
0.2
Corriente B [A]
0
−0.2
−0.4
0.5 1 1.5 0 0.5 1 1.5
7
6 Verdadero
Estimado
5
Posición (Rad)
4
3
2
1
0
0.5 1 1.5 0 0.5 1 1.5
Tiempo [s] Tiempo [s]
Simulación Matriz P
Corriente en el devanado A 0.054 A 0.0936 A
Corriente en el devanado B 0.052 A 0.0936 A
Velocidad 0.26 rad/s 0.4422 rad/s
Posición 0.013 rad 0.0252 rad
Tabla 8.1.: Resultados del ejemplo que muestran una desviación estándar de los
errores de estimación de estado determinados a partir de los resultados
de simulación y determinados a partir de la matriz P del EKF. Estos
resultados son para la simulación del motor de imán permanente de
dos fases. Esta tabla muestra que la matriz P es un buen indicador de
la magnitud de los errores de estimación de estado del EKF.
La matriz P cuantifica la incertidumbre en los estimados de estados. Si las no li-
nealidades en el sistema y las medidas no son demasiado severas, entonces la
matriz P deberá dar una idea de cuan precisa es la estimación. En este ejemplo,
las desviaciones estándar de los errores de estimación de estado se obtuvieron de
la simulación y luego se compararon con los elementos de la diagonal de la matriz
P de estado estacionario que se obtuvo del filtro de Kalman. La Tabla 8.1 muestra
una comparación de los errores de estimación que se determinaron por simulación
y los errores de estimación teóricos basados en la matriz P. Se ve que la matriz
P da una buena indicación de la magnitud de los errores de estimación.
Fk−1 y Lk−1 son definidos por la anterior ecuación. La señal conocida ũk y la
señal de ruido ṽk se definen como:
k , uk , 0) − Fk x̂k
ũk = fk (x̂+ +
(8.28)
ṽk ∼ (0, Lk Qk L⊤
k)
De aquı́,
yk = hk (x̂− −
k , 0) + Hk (xk − x̂k ) + Mk wk
[ ]
= Hk xk + hk (x̂− −
k , 0) − Hk x̂k + Mk wk
= Hk xk + zk + w̃k (8.30)
0 = E {x0 }
x̂+
{ } (8.32)
+ ⊤
P+0 = E (x 0 − x̂+
0 )(x0 − x̂0 )
P− + ⊤ ⊤
k = Fk−1 Pk−1 Fk−1 + Lk−1 Qk−1 Lk−1
(8.34)
x̂− +
k = fk−1 (x̂k−1 , uk−1 , 0)
Kk = P− ⊤ − ⊤
k Hk (Hk Pk Hk + Mk Rk Mk )
⊤ −1
−
[ −
]
k = x̂k + Kk yk − hk (x̂k , 0)
x̂+ (8.36)
P+
k = (I − Kk Hk )P−
k
zk = hk (x̂− −
k , 0) − Hk x̂k
(8.37)
w̃k ∼ (0, Mk Rk M⊤
k)
P− + ⊤ ⊤
k = Fk−1 Pk−1 Fk−1 + Lk−1 Qk−1 Lk−1
Kk = P− ⊤ − ⊤
k Hk (Hk Pk Hk + Mk Rk Mk )
⊤ −1
x̂−
k= fk−1 (x̂+
k−1 , uk−1 , 0)
zk= hk (x̂k , 0) − Hk x̂−
−
k (8.38)
x̂+= x̂k + Kk (yk − Hk x̂−
−
k − zk )
k
[ ]
= x̂k + Kk yk − hk (x̂−
−
k , 0)
−
k = (I − Kk Hk )Pk
P+
ẋ = f (x, u, v, t)
yk = hk (xk , wk )
(8.39)
v(t) ∼ (0, Q)
wk ∼ (0, Rk )
El ruido del proceso v(t) es ruido blanco en tiempo continuo con covarianza Q
8.2. Filtro de Kalman extendido 267
x̂˙ = f (x̂, u, v0 , t)
(8.41)
Ṗ = AP + PA⊤ + LQL⊤
Kk = P− ⊤ − ⊤
k Hk (Hk Pk Hk + Mk Rk Mk )
⊤ −1
−
[ −
]
k = x̂k + Kk yk − hk (x̂k , v0 )
x̂+ (8.42)
P+
k = (I − Kk Hk )P−
k (I
⊤
− Kk H k ) + Kk Mk Rk M⊤ ⊤
k Kk
ẋ = f (x, u, v, t)
yk = hk (xk , wk )
(8.43)
v(t) ∼ (0, Q)
wk ∼ (0, Rk )
0 = E {x0 }
x̂+
{ } (8.44)
+ ⊤
P+0 = E (x 0 − x̂+
0 )(x0 − x̂0 )
x̂˙ = f (x̂, u, 0, t)
(8.45)
Ṗ = AP + PA⊤ + LQL⊤
Kk = P− ⊤ − ⊤
k Hk (Hk Pk Hk + Mk Rk Mk )
⊤ −1
−
[ −
]
k = x̂k + Kk yk − hk (x̂k , 0)
x̂+ (8.46)
P+
k = (I − Kk Hk )P−
k (I
⊤
− Kk Hk ) + Kk Mk Rk M⊤ ⊤
k Kk
y = h(x) (8.48)
define como
[ ]
r
x= (8.49)
θ
Supóngase que x1 (la cual es la magnitud r) es una variable aleatoria con una
media de 1 y una desviación estándar de σr . Supóngase que x2 (el cual es el
ángulo θ) es una variable aleatoria con una media de π/2 y una desviación
estándar de σθ . En otras palabras, las medias de las componentes de x están
dadas como r̄ = 1 y θ̄ = π/2. Además, se asumirá que r y θ son independientes
y que sus funciones de densidad de probabilidad son simétricas alrededor de
sus medias (por ejemplo, gaussiana o uniforme).
Una consideración inicial del anterior problema, junto con la Ecuación (8.47),
podrı́a conducir a creer que y1 tiene una media de 0 y y2 tiene una media de 1.
Además, un enfoque de linealización podrı́a conducir a la misma conclusión. Si
se desarrolla una linealización de primer orden a la Ecuación (8.48) y se toma
el valor esperado de ambos lados, se puede obtener
{ }
∂h
ȳ = E {h(x)} ≈ E h(x̄) + (x − x̄)
∂x x̄
∂h
= h(x̄) + E (x − x̄) = h(x̄)
∂x x̄
[ ]
0
= (8.50)
1
r = r̄ + r̃
(8.51)
θ = θ̄ + θ̃
ȳ1 = 0 (8.53)
ȳ2 = E {r sen θ}
{ }
= E (r̄ + r̃) sen(θ̄ + θ̃)
{ }
= E (r̄ + r̃)(sen θ̄ cos θ̃ + cos θ̄ sen θ̃) (8.54)
No se puede ir más lejos a menos que se asuma alguna distribución para θ̃, ası́
que se asume que θ̃ está uniformemente distribuida entre ±θm . En este caso,
se puede calcular
{ }
ȳ2 = E cos θ̃
sen θm
= (8.56)
θm
Se esperaba obtener 1 por respuesta en confirmación de la Ecuación (8.50), pero
en su lugar se tiene un número que es menor que 1. [Nótese que (sen θm )/θm < 1
para todo θm > 0, y lı́mθm →∞ (sen θm )/θm = 1 ]. El análisis revela un problema
con respecto a la visión inicial y a la linealización de primer orden que se
272 Capı́tulo 8. El filtro de Kalman no lineal
1.01
Media
1 linealizada
0.99
0.98
Media no
2
0.97
y
lineal
0.96
0.95
0.94
0.93
−0.4 −0.3 −0.2 −0.1 0 0.1 0.2 0.3 0.4
y1
realizó antes. La media de y2 será en efecto menor que 1. Esto se puede ver
observando un gráfico de 300 valores de r y t generados aleatoriamente, donde r̃
está uniformemente distribuida entre ±0.01 y θ̃ está uniformemente distribuida
entre ±0.35 rad. La pequeña varianza de r̃ y la gran varianza de θ̃ resulta en
una distribución en forma de arco de puntos como se ve en la Figura 8.2. Esta
distribución en forma de arco resulta en ȳ2 < 1 . Este no es un ejemplo de
filtrado de Kalman. Pero ya que el EKF usa la linealización de primer orden
para actualizar la media del estado, este ejemplo muestra el tipo de error que
puede arrastrarse en el EKF cuando es aplicado a un sistema no lineal.
Para un análisis más general de la media de una transformación no lineal,
recordando de la Ecuación (2.9) que y = h(x) puede ser ampliado en una serie
de Taylor alrededor de x̄ como sigue
1 2 1
y = h(x) = h(x̄) + Dx̃ h + Dx̃ h + Dx̃3 h + · · · (8.57)
2! 3!
8.3. El filtro de Kalman unscented 273
=0 (8.59)
ȳ1 ≈ 0 (8.63)
σθ2 E (θ̃)2
ȳ2 ≈ 1 − =1− (8.64)
2 2
Nótese que se encontró el valor exacto de ȳ2 en la Ecuación (8.55) siendo igual
a E (cos θ̃). La expresión aproximada que se encuentra en la Ecuación (8.64)
corresponde a los dos primeros términos no nulos del desarrollo en serie de
Taylor de E (cos θ̃).
Se puede usar las ecuaciones (8.57) y (8.61) para escribir (y − ȳ) como
[ ] [ ]
1 2 1 1
y − ȳ = h(x̄) + Dx̃ h + Dx̃ h + · · · − h(x̄) + E (Dx̃ h) + E (Dx̃ h) + · · ·
2 4
2! 2! 4!
[ ] [ ]
1 2 1 1
= Dx̃ h + Dx̃ h + · · · − E (Dx̃ h) + E (Dx̃ h) + · · ·
2 4
(8.66)
2! 2! 4!
8.3. El filtro de Kalman unscented 275
sigue:
{ }
⊤ Dx̃ h(Dx̃3 h)⊤ Dx̃2 h (Dx̃2 h)⊤ Dx̃3 h(Dx̃ h)⊤
Py = HPH + E + + +
3! 2! 2! 3!
( 2 ) ( 2 )⊤
Dx̃ h Dx̃ h
+E E + ··· (8.70)
2! 2!
Py ≈ HPx H⊤
[ ][ ][ ] [ 2 ]
0 −1 σr2 0 0 1 σθ 0
= = (8.72)
1 0 0 σθ2 −1 0 0 σr2
1.01
1
Covarianza
0.99 linealizada
0.98
y2
0.97
0.96
0.95 Covarianza
no lineal
0.94
0.93
−0.4 −0.3 −0.2 −0.1 0 0.1 0.2 0.3 0.4
y
1
Esta matriz define una elipse bidimensional, donde Py (1, 1) especı́fica el cua-
drado de la longitud del eje y1 , y Py (2, 2) especı́fica el cuadrado de la longitud
del eje y2 . La Figura 8.3 muestra la covarianza linealizada definida por la
Ecuación (8.72), y la covarianza exacta definida por la Ecuación (8.75). La
covarianza linealizada se centra en la media linealizada y la covarianza exacta
se centra en torno a la media exacta. Se puede observar que la covarianza li-
nealizada no es muy buena aproximación a la covarianza exacta, al menos no
en la dirección de y2 .
Este no es un ejemplo de un filtrado Kalman. pero puesto que el EKF usa
linealización de primer orden para actualizar la covarianza del estado, este
ejemplo muestra la forma del error que puede arrastrarse en el EKF cuando
este es aplicado a un sistema no lineal.
√ (√ )⊤(√ )
donde nP es la raı́z cuadrada de la matriz nP tal que nP nP = nP,
(√ ) √
y nP es la i-ésima raı́z de nP. En las próximas dos subsecciones, se
i
verá como la media del conjunto de los puntos sigma anteriores puede usarse
para aproximar la media y la covarianza de un vector no lineal transformado.
280 Capı́tulo 8. El filtro de Kalman no lineal
Aproximación de la media
Supóngase que se tiene un vector x con una media conocida x̃ y una covarian-
za P, una función no lineal y = h(x), y se desea aproximar la media de y.
Se propone transformar cada punto sigma individual de la Ecuación (8.76)
usando la función no lineal h(·), y entonces se toma la suma ponderada de los
puntos sigma transformados para aproximar la media de y. Los puntos sigma
transformados son calculados como sigue:
( )
y (i) = h x(i) i = 1, ..., 2n (8.77)
∑
2n
ȳu = W (i) y (i) (8.78)
i=1
1 ∑ (i)
2n
ȳu = y (8.80)
2n i=1
Ahora se va a calcular el valor de ȳu para ver lo bien que coincide con la media
real de y. Para hacer esto primero se usa la Ecuación (2.9) para expandir cada
y(i) en la Ecuación (8.80) en una serie de Taylor alrededor de x̄. Esto resulta
en
2n ( )
1 ∑ 1 2
ȳu = h(x̄) + Dx̃(i) h + Dx̃(i) h + · · ·
2n i=1 2!
2n ( )
1 ∑ 1 2
= h(x̄) + Dx̃(i) h + Dx̃(i) h + · · · (8.81)
2n i=1 2!
8.3. El filtro de Kalman unscented 281
=0 (8.82)
Ahora recordando que la media real de y está dada por la Ecuación (8.61)
como
1 { 2 } 1 { 4 }
ȳ = h(x̄) + E Dx̃ h + E Dx̃ h + · · · (8.86)
2! 4!
8.3. El filtro de Kalman unscented 283
1 ∑ 2
∂ h
= E x̃i x̃j
2! i,j=1
∂xi ∂xj x=x̄
1 ∑ ∂ 2 h
n
= E {x̃i x̃j }
2! i,j=1 ∂xi ∂xj x=x̄
1 ∑ ∂ 2 h
n
= Pij (8.87)
2! i,j=1 ∂xi ∂xj x=x̄
Comparando ésta con la Ecuación (8.85) se ve que ȳu (la media aproximada
de y) coincide con la media real de y correctamente hasta el tercer orden,
mientras que la linealización solamente coincide con la media real de y hasta
el primer orden. Si se calcula ȳu usando las ecuaciones (8.76),(8.77) y (8.80),
entonces el valor de ȳu coincidirá con la media real de y hasta el tercer or-
den. La mayor dificultad con este algoritmo es la raı́z cuadrada de la matriz
que es requerida en la Ecuación (8.76). Pero la transformación unscented tiene
la ventaja computacional que la linealización de la matriz H no necesita ser
calculada. Por supuesto, la mayor ventaja de la transformación unscented (re-
lativo a la linealización) es el aumento de la precisión de la transformación de
la media.
Aproximación de la covarianza
∑
2n
Pu = W (i) (y(i) − yu )(y(i) − yu )⊤
i=1
1 ∑ (i)
2n
= (y − yu )(y(i) − yu )⊤ (8.89)
2n i=1
donde los vectores y(i) son los puntos sigma transformados que fueron calcula-
dos en la Ecuación (8.77), y los coeficientes de ponderación W (i) son los mismos
que las dados en la Ecuación (8.79). Expandiendo esta aproximación usando
las ecuaciones (2.9) y (8.83) da lo siguiente:
1 ∑[ ][ ]⊤
2n
Pu = h(x(i) − yu ) h(x(i) − yu )
2n i=1
∑2n [
1 1
= h(x̄) + Dx̃(i) h + Dx̃2 (i) h + Dx̃3 (i) h+ · · ·
i=1
2! 3!
2n ( )]
1 ∑ 1 2 1 4
− h(x̄) − D (j) h + Dx̃(j) h + · · · [· · · ]⊤ (8.90)
2n j=1 2! x̃ 4!
donde
[ ]
1
Υ1 = (Dx̃(i) h) (Dx̃(i) h) = Υ⊤
2 ⊤
1 = 0
2
( )⊤
1 ∑ 1 2
Υ2 = Dx̃(i) h Dx̃(i) h = Υ⊤2 = 0
2n j 2
1 ∑
2n
Pu = (Dx̃(i) h)(Dx̃(i) h)⊤ + TOS (8.92)
2n i=1
donde TOS significa términos de orden superior (es decir, términos a la cuarta
potencia y superior). Expandiendo esta ecuación para Pu dejando de lado los
términos de orden superior da
n ( )( )⊤
1 ∑∑
2n
(i) ∂h(x̄) (i) ∂h(x̄)
Pu = x̃ x̃k (8.93)
2n i=1 j,k=1 j ∂xj ∂xk
cer orden (es decir, solamente los términos a la cuarta y mayor potencia son
incorrectos). Este es el mismo orden de aproximación como el método de linea-
lización. Sin embargo, se espera intuitivamente que la magnitud del error de la
aproximación unscented en la Ecuación (8.89) sea más pequeña que la apro-
ximación lineal HPH⊤ , porque la aproximación unscented al menos contiene
términos indicados correctamente hasta la cuarta potencia y mayor, mientras
que la aproximación lineal no contiene ningún término de otro tipo que HPH⊤ .
La transformación unscented puede resumirse como sigue
La transformación unscented
√ (√ )⊤(√ )
donde nP es la raı́z cuadrada de la matriz nP tal que nP nP =
(√ ) √
nP, y nP es la i-ésima raı́z de nP.
i
3. Transformar los puntos sigma como sigue:
1 ∑ (i)
2n
ȳu = y
2n i=1
1 ∑ (i)
2n
Pu = (y − yu )(y(i) − yu )⊤ (8.97)
2n i=1
8.3. El filtro de Kalman unscented 287
Ejemplo 8.2
Para ilustrar la transformación unscented, considérese la transformación no lineal
mostrada en la Ecuación (8.47). Ya que esas son dos variables independientes
(r y θ), se tiene n = 2. La covarianza de P está dada como P = diag(σr2 , σθ2 ).
La Ecuación (8.79) muestra que W (i)=1/4 para i = 1, 2, 3, 4. La Ecuación (8.76)
muestra que los puntos sigma están determinados como
(√ )⊤
x(1) = x̄ + nP
[ √ ]1
= 1 + σ r 2
π/2
(√ )⊤
x(2) = x̄ + nP
[ 2
]
1 √
=
π/2 + σθ 2
(√ )⊤ (8.98)
x(3) = x̄ − nP
[ √ ]1
= 1 − σr 2
π/2
(√ )⊤
x(4) = x̄ − nP
[ 2
]
1 √
=
π/2 − σθ 2
de y = h(x) como
∑
4
ȳu = W(i) y(i)
i=1
(8.100)
∑4 ( )( )⊤
Pu = W (i)
y(i) − yu y(i) − yu
i=1
xk+1 = f (xk , uk , vk , tk )
(8.112)
yk = h(xk , wk , tk )
8.4. Filtrado de Kalman unscented 289
xk+1 = f (xk , uk , tk ) + vk
yk = h(xk , tk ) + wk
(8.101)
vk ∼ (0, Qk )
wk ∼ (0, Rk )
0 = E {x0 }
x̂+
{ + ⊤
} (8.102)
P0 = E (x0 − x̂+
+
0 )(x0 − x̂0 )
(i)
x̂k−1 = x̂+
k−1 + x̃
(i)
i = 1, . . . , 2n
(√ )⊤
(i) +
x̃ = nPk−1 i = 1, . . . , n
i
(8.103)
(√ )⊤
x̃(n+i) = − nP+
k−1 i = 1, . . . , n
i
Continuación:
(i)
Combinar los x̂k vectores para obtener el estado estimado a priori en
el tiempo k. Esto está basado en la Ecuación (8.80):
1 ∑ (i)
2n
x̂−
k = x̂ (8.104)
2n i=1 k
1 ∑ ( (i) )( )⊤
2n
P− x̂k − x̂− −
(i)
k = k x̂k − x̂k + Qk−1 (8.105)
2n i=1
x̂k = x̂−
(i) (i)
k + x̃ i = 1, . . . , 2n
(√ )⊤
(i) −
x̃ = nPk i = 1, . . . , n
i
(8.106)
(√ )⊤
x̃(n+i) = − nP−
k i = 1, . . . , n
i
Si se desea, este paso puede omitirse. Esto es, en vez de generar nue-
vos puntos sigma se pueden reutilizar los puntos sigma que fueron
obtenidos a partir de la actualización de tiempo. Esto ahorrará es-
fuerzo computacional si se está dispuesto a sacrificar el rendimiento.
Usar la ecuación no lineal de medición h(·) para transformar los
(i)
puntos sigma en ỹk vectores (medida prevista), como se muestra
en la Ecuación (8.76):
(i) (i)
ŷk = h(x̂k , tk ) (8.107)
8.4. Filtrado de Kalman unscented 291
Continuación:
(i)
Combinar los ŷk vectores para obtener la medida prevista en el tiempo
k. Esto está basado en la Ecuación (8.80):
1 ∑ (i)
2n
ŷk = ŷ (8.108)
2n i=1 k
1 ∑ ( (i) )( )⊤
2n
(i)
Py = ŷk − ŷk ŷk − ŷk + Rk (8.109)
2n i=1
1 ∑ ( (i) )( )⊤
2n
x̂k − x̂−
(i)
Pxy = k ŷ k − ŷk (8.110)
2n i=1
Kk = Pxy P−1
y
−
k = x̂k + Kk (yk − ŷk )
x̂+ (8.111)
P+
k = P−
k − Kk P y K⊤
k
como en []:
xk
= vk
(a)
xk (8.113)
wk
(a)
Entonces se puede utilizar el UKF para estimar el estado aumentado xk . El
UKF se inicializa como
E {x0 }
x̂0 = 0
(a)
0
{ } (8.114)
E (x0 − x̂0 )(x0 − x̂0 )⊤ 0 0
P0 = Q0 0
(a)+
0
0 0 R0
ρ0 = 2 lb − s/f t4
g = 32.2 f t/s2
k = 20.000 f t
{ 2}
E wk = 10.000 f t2 (8.116)
{ }
E vi2 (t) = 0 i = 1, 2, 3
M = 100.000 f t
a = 100.000 f t
Las condiciones iniciales del sistema y el estimador están dadas como
⊤
x0 = [300.000 −20.000 0.001]
x̂+
0 = x0
[ ] (8.117)
1.000.000 0 0
+
P0 = 0 4.000.000 0
0 0 10
Se ha visto que una media precisa y una covarianza aproximada para una
transformación no lineal y = h(x) se puede obtener escogiendo 2n puntos
sigma (donde n es la dimensión de x) como se da en la ecuación ,la media y
la covarianza aproximadas como se da en las ecuaciones y . Sin embargo, se
puede demostrar que el mismo orden de precisión de la estimación de la media
y la covarianza se puede conseguir mediante la elección (2n + 1) puntos sigma
x(i) como sigue
x(0) = x̄
x(i) = x̄ + x̃(i) i = 1, . . . , 2n
(√ )⊤
(8.118)
x̃(i) = (n + κ)P i = 1, . . . , n
i
(√ )⊤
x̃(n+i) = − (n + κ)P i = 1, . . . , n
i
1. Escoger el peso W (0) ∈ [0, 1). La escogencia de W (0) afecta sólo a los
momentos de cuarto y de orden superior del conjunto de puntos sigma[].
2W (1) (8.122)
1
σ2 = √
(1)
2W (1)
1. Escoger el peso W (0) ∈ [0, 1). La escogencia de W (0) afecta sólo a los
momentos de cuarto y de orden superior del conjunto de puntos sigma[].
1 − W (0)
Wi = i = 1, . . . , n + 1 (8.125)
n−1
Note que (en contraste a la transformación unscented sı́mplex) todos los
pesos son idénticos excepto para W (0) .
2W (1) (8.126)
1
σ2 = √
(1)
2W (1)
(n)
La relación del elemento más grande de σi al elemento más pequeño es
/
n 1
√ √ =n (8.129)
n(n + 1)W (1) n(n + 1)W (1)
ası́ los problemas numéricos no podrán ser tema para la transformación uns-
cented esférica.
Example : Considérese el sistema de un cuerpo cayendo como el descrito en el
ejemplo[]. Las condiciones iniciales del sistema y el estimador están dadas como
[ ]⊤
x0 = 300.000 −20.000 1/1000
[ ]⊤
0 = 303.000 −20.200 1/1010
x̂+
(8.130)
30.000 0 0
P+0 =
0 2.000 0
0 0 1/10.000
Se corren 100 simulaciones Monte Carlo, cada una con un tiempo de simulación
de 60 s. El promedio de errores de estimación RMS del EKF, el estándar UKF
(seis puntos sigma), el UKF sı́mplex (cuatro puntos sigma puesto que se escogió
W (0) = 0), y el UKF espiral (cuatro puntos sigma puesto que se escogió W (0) = 0)
están dados en la tabla[]. El UKF sı́mplex funciona mejor para estimación de
altitud, con el UKF estándar no muy lejos. El UKF estándar se desempeña mejor
para la estimación de velocidad, y el UKF espiral se desempeña mejor para la
estimación del coeficiente balı́stico. El EKF es en general el de peor desempeño
de los cuatro estimadores de estado.
300 Capı́tulo 8. El filtro de Kalman no lineal
Athans, M. and Safonov, M., G. (1978). Robustness and computational aspects of nonlinear
stochastic estimators and regulators. IEEE Transactions on Automatic Control, AC-
23(4): 717–725.
Bellantoni, J. F. and Dodge, K. W. (1967). A square root formulation of the Kalman–
Schmidt filter. AIAA Journal, 5(7): 1309–1314.
IEEE (2012). First-hand: The unscented transform. Global History Network. URL:
http://www.ieeeghn.org/wiki/index.php/First-Hand:The_Unscented_Transform.
Simon, D. (2006). Optimal State Estimation: Kalman, H-infinity, and Nonlinear Approaches.
John Wiley & Sons, Inc., New Jersey, USA.
301
Apéndice A
Funciones matriciales
∂Tr(F(X))
= f ⊤ (X)
∂X
donde f (·) es la derivada escalar de F(·).
Primer Orden
∂
Tr(X) = I (A.1)
∂X
∂
Tr(XA) = A⊤ (A.2)
∂X
∂
Tr(X⊤ A) = A (A.3)
∂X
∂
Tr(AX⊤ ) = A (A.4)
∂X
∂
Tr(AXB) = A⊤ B⊤ (A.5)
∂X
303
304 Capı́tulo A. Funciones matriciales
∂
Tr(AX⊤ B) = BA (A.6)
∂X
∂
Tr(A ⊗ X) = Tr(A)I (A.7)
∂X
Segundo Orden
∂
Tr(X2 ) = 2X⊤ (A.8)
∂X
∂
Tr(X2 B) = (XB + BX)⊤ (A.9)
∂X
∂ ∂
Tr(X⊤ X) = Tr(XX⊤ ) = 2X (A.10)
∂X ∂X
∂
Tr(X BX) = BX + B⊤ X
⊤
(A.11)
∂X
∂
Tr(BXX⊤ ) = BX + B⊤ X (A.12)
∂X
∂
Tr(XX⊤ B) = BX + B⊤ X (A.13)
∂X
∂
Tr(XBX⊤ ) = XB⊤ + XB (A.14)
∂X
∂
Tr(BX⊤ X) = XB⊤ + XB (A.15)
∂X
∂
Tr(X⊤ XB) = XB⊤ + XB (A.16)
∂X
∂
Tr(AXBX) = A⊤ X⊤ B⊤ + B⊤ X⊤ A⊤ (A.17)
∂X
∂
Tr(B⊤ X⊤ CXB) = C⊤ XBB⊤ + CXBB⊤ (A.18)
∂X
∂
Tr(X⊤ BXC) = BXC + B⊤ XC⊤ (A.19)
∂X
∂ [ ]
Tr (AXB + C)(AXB + C)⊤ = 2A⊤ (AXB + C)B⊤ (A.20)
∂X
∂ ∂
Tr(X ⊗ X) = Tr(X)Tr(X) = 2Tr(X)I (A.21)
∂X ∂X
A.1. Derivada de la traza 305
Orden Superior
∂
Tr(Xk ) = k(Xk−1 )⊤ (A.22)
∂X
∂ ∑
k−1
Tr(AXk ) = (Xr AXk−r−1 )⊤ (A.23)
∂X r=0
∂ [ ⊤ ⊤ ]
Tr B X CXX⊤ CXB = CXX⊤ CXBB⊤ + C⊤ XBB⊤ X⊤ C⊤ X
∂X
+ CXBB⊤ X⊤ CX + C⊤ XX⊤ C⊤ XBB⊤
(A.24)
Otras
∂
Tr(AX−1 B) = −(X−1 BAX−1 )⊤ = −X−⊤ A⊤ B⊤ X−⊤ (A.25)
∂X
Supóngase que B y C son simétricas, entonces
∂ [ ]
Tr (X⊤ CX)−1 A = −(CX(X⊤ CX)−1 )(A + A⊤ )(X⊤ CX)−1
∂X
(A.26)
∂ [ ]
Tr (X⊤ CX)−1 X⊤ BX = −2CX(X⊤ CX)−1 X⊤ BX(X⊤ CX)−1
∂X
−1
+ 2BX(X⊤ CX) (A.27)
∂ [ ]
Tr (A + X⊤ CX)−1 X⊤ BX = −2CX(A + X⊤ CX)−1 X⊤ BX(A + X⊤ CX)−1
∂X
−1
+ 2BX(A + X⊤ CX) (A.28)
Apéndice B
Transformada delta
qxk = Aq xk + Bq uk (B.2)
yk = C q xk (B.3)
∆t q−1
δ= (B.4)
T
307
308 Capı́tulo B. Transformada delta
δxk = Aδ xk + Bδ uk , (B.6)
yk = Cδ xk . (B.7)
q = 1 + Tδ. (B.8)
donde
∫ T
1 AT A2 T2
Λ= eAτ dτ = I + + + ··· (B.11)
T 0 2! 3!
Este resultado se genera directamente a partir de la consideración de una
expansión en series de potencia de Aδ = (eAT − I)/T. Para determinar Λ
se aplican métodos tradicionales, tales como autovalores, teorema de Cayley–
Hamilton, etc. A partir de (B.10) y (B.11), se ve que hay una relación cercana
entre los modelos en el dominio delta y los modelos en el dominio del tiempo
continuo, dado que Λ → I cuando ∆t → 0.
donde
vq
vδ = , wδ = vq (B.14)
T
De aquı́ se obtiene para las covarianzas
Qq
cov{vδ } = , Qq (B.15)
T2
cov{wδ } = Rq , Rq (B.16)
Sq
cov{vδ , wδ } = , S′ δ ; Sδ = TS′ δ (B.17)
T
Para muestreo rápido Qq , Rq son tı́picamente de orden T y 1/T, respectiva-
mente. Ası́, Qδ , Rδ son ambos del orden 1/T. Por lo tanto, para mantener
310 Capı́tulo B. Transformada delta
Qδ Qq
Ωδ , = TQδ = (B.18)
fs T
Rδ
Γδ , = TRδ = TRq (B.19)
fs
Estas relaciones pueden ser utilizadas para desarrollar las formas con operador
delta de los algoritmos del filtro de Kalman desde sus formas con funciones de
desplazamiento.
⊤ −1
= P+k Hk Rk
− −
k = x̂k + Kk (yk − Hk x̂k ) = estado estimado a posteriori
x̂+
− ⊤ ⊤
(B.20)
k = (I − Kk Hk )Pk (I − Kk Hk ) + Kk Rk Kk
P+
= [(P−k)
−1
+ H⊤ −1
k Rk Hk ]
−1
= (I − Kk Hk )P−k
Las ecuaciones anteriores se pueden separar en dos grupos: las primeras dos,
predicen el estimado del estado y la matriz de covarianza del error desde los
valores anteriores y pueden, como tales, ser denominadas como ecuaciones de
predicción. Estas predicciones se usan para formar la ganancia de Kalman, la
cual es usada en las últimas lı́neas junto con la última muestra de salida en la
forma de una señal de innovación, para corregir las predicciones y producir un
estimado a posteriori de xk y Pk . Éstas pueden entonces, denominarse como
las ecuaciones de corrección.
Tomando la relación entre el operador delta y el operador de desplazamiento
en las descripciones de espacio de estado mostradas en (B.12) a (B.17), y
entonces sustituyéndolas en las ecuaciones (B.20), el filtro de Kalman con el
operador delta se puede construir en su forma de predicción corrección de la
B.1. Forma delta del estimador de estado óptimo 311
siguiente manera:
x̂−
k = (T δFk + I)x̂− δ
k−1 + T Gk uk−1 (B.21)
δ −
Pk = (T δFk + I)x̂− δ ⊤ 2δ
k−1 (T Fk + I) + T Qk (B.22)
δ −δ ⊤ δ − ⊤
δ
Kk = Pk Hk [ Hk δPk δHk + δRk ]−1 (B.23)
x̂+
k = x̂− δ δ −
k + Kk [yk − Hk x̂k ] (B.24)
+ −
δ
Pk = (I − δKk δHk ) δPk (B.25)
313