Вычислительные методы алгебры и оценивания. И.В. Семушкин (2011) (1185350), страница 43
Текст из файла (страница 43)
Обработка измерения:e=xb(t−А. Начальное присваивание: Pe = P (t−i ).i ), xБ. m-кратное повторение процедуры скалярного обновления.Для j = 1, 2, . . . , m выполнять:α := hT Peh + rj (ti ) , v := Peh , K := v/α ,Pb := Pe − Kv T , v := Pbh , Pb := Pb − vK T + KK T ,xb := xe + K(z − hT xe)с экстраполяцией между повторениями: Pe := Pb ,xe := xb.27513 Устойчивые алгоритмы фильтрацииbb,xb(t+В. Завершающее присваивание: P (t+i ) := xi ) := P ,где h — j-й столбец матрицы H T (ti) z — j-й элемент вектора z(ti ) ,j = 1, 2, .
. . , m .13.5Квадратно-корневой фильтр ПоттераКак отмечалось в подразд. 11.8, здесь вместо матриц P (t±i ), по своей природе положительно определенных, оперируют с их квадратными корнями±±T ±S(t±i ), отвечающими равенствам S(ti )S (ti ) = P (ti ). Ясно, что эти соотношения не дают однозначного определения квадратных корней. Однакоэто обстоятельство в общем случае не вызывает беспокойства, посколькуэти корни однозначно определяются по разложению Холесского (например,LLT -разложение или другие, см. разд.
6).Основная идея метода фильтрации с использованием квадратного корнясостоит в замене уравнений алгоритма Калмана на аналогичные, предназначенные для последовательного расчета S(t±i ). Такой подход оправдывается±T ±тем, что произведение S(ti )S (ti ) не теряет свойство положительной определенности (при условии полноты ранга) даже с учетом ошибок округления,тогда как ошибки округления могут приводить к потере этого свойства дляматрицы P (t+ ), если она вычисляется по стандартному алгоритму (13.3).На этапе экстраполяции, с учетом разложения матрицы P (t−i ) и также1/2T /2Q(ti−1) = Q (ti−1)Q (ti−1), уравнение для нее принимает следующий вид:T −+T +TS(t−i )S (ti ) = Φ(ti , ti−1 )S(ti−1)S (ti−1)Φ (ti , ti−1 )++ Γ(ti−1)Q1/2(ti−1)QT /2(ti−1)ΓT (ti−1).Непосредственно вычислять матрицу S(t−i ) можно путем построения ортогональной матрицы T размера (n + s) × (n + s), такой что T + T − S (ti−1)ΦT (ti, ti−1)S (ti )=T.0QT /2(ti−1)ΓT (ti−1)Это обеспечивает любой алгоритм ортогональных преобразований.
Хорошийрезультат дает [80] модифицированный алгоритм Грама–Шмидта. В результате после этапа экстраполяции матрица S(t−i ) всегда будет получаться треугольной.На этапе обработки измерения требуем иметь уравнение Pb := Pe − KhPee n − f f T /α)SeT . Определим число β так, чтобы обеспечитьв виде SbSbT := S(Iтождественность: In − f f T /α = (In − βf f T )(In − βf f T ).27613.5 Квадратно-корневой фильтр ПоттераИз получающегося отсюда квадратного уравнения, с учетом диагональности R(ti ) = diag [r1(ti ), r2(ti), .
. . , rm(ti )] и равенства α = f T f + rj (ti ), выбираем одно решениеpβ = (1/α)/(1 + 1/α),pзащищенное от операции вычитания положительной величины 1/α в знаменателе.С учетом вышеизложенного получим алгоритм фильтра Поттера:I. Экстраполяция:xb(t−x(t+i ) = Φ(ti , ti+1 )bi−1 );S T (t−i )0=Txb(t+0 ) = x̄0 ;TS T (t+i−1)Φ (ti , ti−1 )QT /2(ti−1)ΓT (ti−1).II. Обработка измерения:e=xb( t−А. Начальное присваивание: Se = S(t−i ).i ); xБ. m-кратное повторение процедуры скалярного обновления.Для j = 1, 2, . .
. , m выполнять:f := SeT h;α := f T f + rj (ti ); γ := 1/(1 +e /α; Sb := Se − γKf T ;K := Sfxb := xe + K(z − hT xe)bс экстраполяцией между повторениями Se := S;p1/α);xe := xb.b xВ. Завершающее присваивание: S(t+b(t+b,i ) := S;i ) := xTгде h — j-й столбец матрицы H (ti ); z — j-й элемент вектора z(ti ) ,j = 1, 2, . . .
, m.Теперь вычисление S равносильно счету с двойной точностью P в стандартном алгоритме Калмана. Кроме того, устранена опасность утраты матрицей Pb = SbSbT свойства положительной определенности, что влекло бырасходимость оценок вектора состояния. Недостатком данного алгоритмаявляется для каждогоp момента ti m-кратное наличие операции извлеченияквадратного корня 1/α.27713 Устойчивые алгоритмы фильтрации13.6Одноранговое обновление ПО-матрицРассмотрим теорему об одноранговом обновлении положительно определенной матрицы [15], которую доказали Agee и Turner (1972) для версииU DU T разложения ПО–матрицы.
Переформулируем ее для варианта LDLTразложения.Теорема 13.1 (Agee–Turner, 1972). Пусть матрица L — нижняя треe = diag [de1, de2, . . . , den ] — диаугольная с единичной диагональю, матрица De T . Пусть также заданы некоторый скаляр c и вектор–гональная и P = LDLстолбец a = [a1 , a2, . . . , an ]T , такие что 0 < P̄ = P + caaT . Тогда разложениеP̄ = L̄D̄L̄T , где L̄ — нижняя треугольная с единичной диагональю матрицаи D̄ = diag [d¯1, d¯2, . . .
, d¯n], существует и дается следующим алгоритмом:А. Начальное присваивание: c1 = c.Б. Для i = 1, 2, . . . , n − 1 выполнять:1) d¯i = dei + ci a2i ;2) Для k = i + 1, i + 2, . . . , n выполнять:i. ak := ak − ai lki ;ii. ¯lki = lki + ci ai ak /d¯i; ✮ В матрицах L нетривиальные элементынаходятся только ниже диагонали.3) ci+1 = ci dei /d¯i.В. Завершающее присваивание: d¯n = den + cn a2n .Доказательство.
Запишем квадратичную форму 0 < xT P̄ x в видесуммы полных квадратов с подстановкой в нее левой и правой частей равенства P̄ = P + caaT . Продолжим доказательство в виде несложного упражнения, последовательно выделяя эти полные квадраты с правой стороныполученного равенства квадратичных форм (см.
[78]). Проделайте это самостоятельно в виде упражнения.2Упражнение 13.1.версию Теоремы 13.1.Аналогично предыдущему, докажите следующуюe — нижняя треугольная матрица и Pe = LeLeT .Теорема 13.2. Пусть LСкаляр c и вектор–столбец a = [a1 , a2 , . . . , an ]T такие что 0 < P̄ = Pe + caaT .Тогда разложение P̄ = L̄L̄T существует и дается следующим алгоритмом:27813.7 Факторизованный фильтр БирманаА.
Начальное присваивание: c1 = c .Б. Для i = 1, 2, . . . , n выполнять:q¯1) lii = el 2ii + ci a2i ;2) Для k = i + 1, i + 2, . . . , n выполнять:eei. ak := ak − ai lki/lii ;ee¯ii. lki = lki /lii ¯lii + ci ai /¯lii ak ;2e3) ci+1 = ci l ii /¯l 2ii.Упражнение 13.2.e=L2 0−1 3,a=1−22,c = −1 .√3030eLeT + caaT =√Найдите P̄ = L, L̄ =. Алгоритм тео0 606ремы 13.2 дает этот же результат.
Проверьте. Измените исходные данные, например, возьмитечто алгоритм дает правильный√ c = 1. Проверьте,p 0 .√5результат: L̄ =−4/ 53 6/5Упражнение 13.3. Сформулируйте и докажите еще две версии теоремыоб одноранговом обновлении подобно двум предыдущим теоремам, опираясьна другие разложения Холесского: P = U U T и P = U DU T [15].13.7Факторизованный фильтр БирманаОсновная идея этого алгоритма состоит в разложении ковариационной матрицы P в произведение двух треугольных матриц и диагональной матрицы между ними.
Можно рассматривать два варианта алгоритма: LD-алгоритм или U D-алгоритм Бирмана. Последний изложен в[15]. Здесь рассмотрим LD-алгоритм, т. е. используем разложение P (t±i ) =±±±T ±= L(ti )D(ti )L (ti ), где L(ti ) — нижние треугольные матрицы с единичнойдиагональю, D(t±i ) — диагональные матрицы.27913 Устойчивые алгоритмы фильтрацииНа этапе экстраполяции с учетом разложения матрицы Q(ti−1) == Lq (ti−1)Dq (ti−1)LTq (ti−1) и аналогичного разложения матрицы P (t−i ) урав−нение для P (ti ) принимает следующий вид:−−T −++T +TP (t−i ) = L(ti )D(ti )L (ti ) = Φ(ti , ti−1 )L(ti−1)D(ti−1)L (ti−1)Φ (ti , ti−1 ) ++ Γ(ti−1)Lq (ti−1)Dq (ti−1)LTq (ti−1)ΓT (ti−1).Представим матрицу P (t−i ) в следующем виде:+ Γ(ti−1)Lq (ti−1) ×))=Φ(t,t)L(tP (t−ii−1ii−1 T +T)Φ(t,t)L(ti i−1i−1× Diag [D(t+,TTi−1), Dq (ti−1)]Lq (ti−1)Γ (ti−1)Tт.
е. P (t−i ) = W (ti )D(ti )W (ti ),Tw(t)i1 Γ(ti−1)Lq (ti−1) = . . . ,)W (ti) = Φ(ti , ti−1)L(t+i−1wnT (ti )D(ti ) = Diag [D(t+i−1), Dq (ti−1)] = Diag [D1 (ti ), . . . , DN (ti )],где W (ti) – матрица размера (n×(n+s)) и N = n+s, s — размерность шума−w(ti ) в уравнении состояния (13.1). Факторы L(t−i−1) и D(ti−1 ) вычисляютсяметодом модифицированного взвешенного алгоритма Грама–Шмидта [15] наэтапе экстраполяции.Перефразируя его для нижних треугольных факторов, запишем следуe L,b D,e D,b D, Wющий результат (для удобства изложения будем писать L,−+−+вместо L(ti ), L(ti ), D(ti ), D(ti ), D(ti ), W (ti) соответственно):Теорема 13.3 (Взвешенная ортогонализация Грама–Шмидта и факторизация матрицы — аналог теоремы VI.4.1 из [15]).Пусть векторыNw1, .
. . , wn образуют линейно независимую систему, wi ∈ R , N ≥ n и D == Diag [D1, . . . , Dn] > 0. Применим следующий алгоритм (нижний индексуказывает номер нетривиального элемента матрицы / вектора):(1)А. Начальное присваивание: vk = wk , k = 1, . . . , n.Б. Для j = 2, . .
. , n выполнять:280(j−1)(j−1)dej−1 = (vj−1 )T Dvj−1 ;(j−1)(j−1)elk,j−1 = (vk )T Dvj−1 /dej−1, k = j, . . . , n;(j−1)(j−1)(j)−elk,j−1vj−1 , k = j, . . . , n .vk = vk13.7 Факторизованный фильтр Бирмана(n)(n)В. Завершающее присваивание: den = (vn )T Dvn .(j)Тогда vj = vj , j = 1, . . . , n, где vj — взвешенные взаимно ортогональныевекторы: viT Dvj = dej δi,j (δi,j — символ Кронекера) и T T v1w1 Te =LeDeLeT .e . .
. D v1 . . . vn L . . . D w1 . . . wn = LvnTwnT−eeИз этого алгоритма имеем результат: L(t−i ) = L, D(ti ) = D.Этап обработки измерений представим в виде теоремы, формулируя еедля нижних треугольных сомножителей (LD-версия) [78] вместо U D-версиитеоремы Бирмана из [15].Теорема 13.4 (Алгоритм Бирмана). Пусть калмановская процедураbDbLbT искалярного обновления (13.5), (13.6) использует разложения Pb = LeDeLeT , где L — нижние треугольные с единичной диагональю матрицы,Pe = LD — диагональные (с положительными элементами) матрицы. Тогда даннаяпроцедура (13.5), (13.6) эквивалентна пунктам В, Г, Д и Е, следующегоалгоритма.II. Обработка измерения:e = D(t− ), xe = L(t− ), De = x(t−А. Начальное присваивание: Li ).iiБ. m-кратное повторение процедуры скалярного обновления.Для j = 1, 2, .
. . , m выполнять:eT h;В. Вычислить векторы f = [f1, f2, . . . , fn]T = Le .v = [v1, v2, . . . , vn]T = DfГ. Задать начальные значения α 0 = r; K = [0 . . . 0 vn ]T .Д. Для i = n, n − 1, . . . , 2, 1 выполнять:началоα := α 0 + vifi ; γ := 1/α ;dbi := deiα 0 γ ;λ := −fi γ ;bli := eli + λK ; K := K + eli v i ;(✮)0α := α .конецЕ. Вычислить векторы ν := γ(z − hT xe) ; xb := xe + Kνс экстраполяцией между повторениями:e := Lb; De := Db; xLe := xb.28113 Устойчивые алгоритмы фильтрацииЖ.
Завершающее присваивание:bbL(t+D(t+xb(t+b.i ) := L ;i ) := D ;i ) := xЗдесь h — j-й столбец матрицы H T (ti ); z — j-й элемент вектораz(ti ); r — j-й элемент rj (ti ) диагональной матрицы ковариацийшума измерений R(ti), j = 1, 2, . . . , m — номер скалярного измерения в составе вектора измерений z(ti ) в момент ti .Замечание 13.5. Строка ✮ в алгоритме на стр. 281 выглядит так:Для k = i + 1 до n выполнятьначалоblki := elki + λKk , Kk := Kk + elkiviконецОна пропускается при i = n. Здесь Kk есть k-й элемент того вектора K,который существует в цикле Д алгоритма на стр. 281.Доказательство.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
