Вычислительные методы алгебры и оценивания. И.В. Семушкин (2011) (1185350), страница 44
Текст из файла (страница 44)
Применим для (13.5) факторизацию Холесского видаbDb PbT и Pe = LeDeLe T . ПолучимPb = Lс обозначениямиbDb PbT = L(e De + cvv T )LeTLe T h, v = Df.ef =Le + cvv T , в ней следует считать, что:Применяя теорему 13.1 к выражению De T + cvv T приP ≡ D, c ≡ −1/α и a ≡ v. Отсюда получаем L̄D̄ L̄T = LDLтом, что в этом выражении имеем L = I. Когда L̄ и D̄ будут найдены, мыполучимbDb PbT = LeL̄D̄ L̄T LeT =⇒ Db = D̄, Lb=LeL̄.L(13.7)e + cvv T . Получим:Запишем теорему 13.1 применительно к выражению Dc = −1/α,А. Начальное присваивание: c1 = c.Б. Для i = 1, 2, .
. . , n − 1 выполнять:1) dbi = dei + ci vi2. ①2) Для k = i + 1, i + 2, . . . , n выполнять:282i. vk := vk − vilki = vk . ❍ Здесь L = I, т. е. lki = 0. Это действиеисчезает.ii. ¯lki = lki + (civi /dbi)vk = (ci vi/dbi)vk . ②13.7 Факторизованный фильтр Бирмана3) ci+1 = ci dei /dbi. ③В. Завершающее присваивание: dbn = den + cn vn2 .④Этот алгоритм теоретически верен, но численно неустойчив из-за того,что должно выполняться условие: ∀i : ci < 0. В силу этого величины dbiв ① могут стать отрицательны, что противоречит требованию положительной определенности Pb > 0. Устраним этот недостаток, эквивалентно преобразуя алгоритм на стр.
282–283 по пунктам ①, ②, ③, ④. Преобразуем ①:dbi /ci = dei /ci + vi2. С учетом ③: dei/ci+1 = dbi/ci , поэтому dei/ci+1 = dei/ci + vi2 ,следовательно, 1/ci+1 = 1/ci + vi2 /dei. Однако vi = dei fi. Поэтому −1/ci == −1/ci+1 + vi fi. Введем обозначение: ∀i : αi = −1/ci. Это позволяет записать:αi = αi+1 + vifi , i = n − 1, n − 2, . . .
, 2, 1.(13.8)Финальное значение, т. е. α1 , согласно введенному обозначению, должно совпасть с α = r + v T f , где r = r(tj ), так как c1 = c и c = −1/α. Чтобы этопроизошло, в алгоритме необходимо стартовать от значения αn = r + vnfn .Таким образом, из ① получили алгоритм (13.8), заменяющий ③.Теперь выведем из ③ алгоритм, заменяющий ①:dbi = dei (ci /ci+1) = dei (αi+1/αi ) , i = n − 1, n − 2, . . . , 2, 1.Осталось преобразовать ②. Из ② имеем ¯lki = λi vk , λi = (ci vi/dbi), а из③ ci /dbi = ci+1/dei, то есть λi = −fi/αi+1, так как vi = dbi fi. Следовательно,алгоритм для ¯lki вместо ② приобретает вид:Для i = n − 1, n − 2, .
. . , 2, 1 вычислятьλi = −fi /αi+1 .Для k = i + 1 до n вычислять¯lki = λi vk .Наконец, преобразуем последнее действие ④ на стр. 283. Последовательнонаходим:dbn /cn = den /cn + vn2 ,(dbn /den)(1/cn) = (1/cn) + (vn2 /den ),(v 2 /den ) = (den fn)2 /den = vn fn .nИз обозначения (−1/ci+1) = αi+1 при i = n − 1 следует (−1/cn) = αn , т. е.αn = −(dbn /den)(1/cn) + vn fn. С другой стороны, известно стартовое значение αn = r + vnfn , следовательно, dbn = den r/αn . Запишем результирующийэквивалентный алгоритм:28313 Устойчивые алгоритмы фильтрацииЗадать начальные значения: αn = r + vn fn, dbn = denr/αn .Для i = n − 1, n − 2, .
. . , 2, 1 вычислятьαi = αi+1 + vi fi ;dbi = dei (αi+1/αi ) ;λi = −fi/αi+1 .Для k = i + 1 до n вычислять¯lki = λi vk .Отмеченный выше недостаток устранен: теперь ∀i : dbi > 0. Проведенное эквивалентное преобразование привело к инверсии направления вычисe и Deлений. В исходном массиве для нетривиальных элементов матриц Lисходные столбцы сначала имеют (начиная с диагонали вниз) вид:"#"#"#hieeed1d2dn−1e,, ...
, e,dn .eel1l2ln−1После выполнения этого алгоритма они имеют следующее наполнение:hidbn−1db2db1bdn .¯l2 , . . . , ¯ln−1 ,¯l1 ,Вернемся к последней формуле в (13.7), чтобы в согласии с ней перейтиbот «промежуточной» матрицы L̄ к «окончательной» матрице L.Матрица L̄ — нижняя треугольная с единицами на главной диагонали.В ней поддиагональная часть i-го столбца, согласно нижней строке алгоритма на стр. 284, найдена в виде ¯liT = λi [vi+1, vi+2, . .
. , vn], где использованы элементы из вектора v T = [v1, v2, . . . , vn], начиная от элемента vi+1 идалее. Введем вспомогательные обозначения n-мерных векторов:000v1 .. .. v v . . 2 2 (n−1) v (1) v (n)(2)v = 0 , v= 0 , . . . , v = 3 , v = 3 = v. .. .. 0 vn−1 . . vnvnvnvnТеперь видно, чтоhL̄ = I + λ1 v284(2) i λ2 v (3) · · · λn−1 v (n) 0 .13.8 Квадратно-корневой фильтр КарлсонаОтсюдаh i(3) (2) (n) ebeeee· · · λn−1Lvλ2 LvL = LL̄ = L + λ1 Lv0 .Введем обозначения возникающих здесь матриц-столбцов:e (2) ,K2 = Lve (3) , · · · , Kn = Lve (n) .K3 = Lve (n) .
ДалееВ качестве начального столбца имеем Kn = Lve (n−1) = Kn + eKn−1 = Lvln−1vn−1 ,e (n−2) = Kn−1 + eKn−2 = Lvln−2vn−2 ,··· ,Ki = Ki+1 + eli vi, i = n − 1, n − 2, . . . , 2, 1 ,и в конце получаемe (1) = Lve = K̄,K1 = LvK̄ — вспомогательное обозначение, причем здесь для строгости записей надопонимать, что eli , i = n−1, n−2, . .
. , 2, 1, обозначает весь j-й столбец матрицыe — вместе с тривиальными элементами 0 и 1.L2Данный в теореме 13.4 алгоритм не содержит операции извлечения квадратного корня, а работа с треугольными матрицами требует меньшего числаарифметических операций по сравнению с обычными.13.8Квадратно-корневой фильтр КарлсонаЭтап экстраполяции здесь в точности совпадает с этим этапом в алгоритме Поттера. Выведем только алгоритм этапа обработки измерения.Упражнение 13.4. Выведите алгоритм Карлсона обработки измерения. Вывод можно сделать в полном соответствии с тем, как выше быладоказана теорема 13.4 (см. стр.
281). Другой способ вывода — как следствиеe принять формально то, что там было LeDe 1/2. Первый спотеоремы 13.4: за Lсоб предпочтительнее для понимания доказательства, второй полезен дляосвоения техники вычислений. Получите следующий результат.Теорема 13.5. Пусть калмановская процедура скалярного обновлеbLbT и Pe = LeLeT , где Lb иния (13.5), (13.6) использует разложения Pb = LL — нижние треугольные матрицы. Тогда данная процедура (13.5), (13.6)эквивалентна пунктам В, Г, Д и Е, следующего алгоритма.28513 Устойчивые алгоритмы фильтрацииII. Обработка измерения:e = L(t− ), xe = x(t−А. Начальное присваивание: Li ).iБ.
m-кратное повторение процедуры скалярного обновления.Для j = 1, 2, . . . , m выполнять:eT h .В. Вычислить векторы f = [f1, f2, . . . , fn]T = LГ. Задать начальные значения: α 0 = r; K = [0 . . . 0 elnn fn ]T .Д. Для i = n, n − 1, . . . , 2, 1 выполнять:началоpα := α 0 + fi2 ;β :=α 0 /α ;blii := βelii ;λ := −fi/(αβ) ;bli := βeli + λK ; K := K + eli f i ;(✮)α 0 := α .конецЕ.
Вычислить векторы ν := (z − hT xe)/α ; xb := xe + Kνс экстраполяцией между повторениями:e := Lb; xLe := xb.Ж. Завершающее присваивание:bb.xb(t+L(t+i ) := xi ) := L ;Здесь h — j-й столбец матрицы H T (ti ); z — j-й элемент вектораz(ti ); r — j-й элемент rj (ti ) диагональной матрицы ковариацийшума измерений R(ti), j = 1, 2, . . . , m — номер скалярного измерения в составе вектора измерений z(ti ) в момент ti .Замечание 13.6.
Строка ✮ в алгоритме на стр. 286 выглядит так:Для k = i + 1 до n выполнятьначалоblki := βelki + λKk , Kk := Kk + elkifiконецОна пропускается при i = n. Здесь Kk есть k-й элемент того вектора K,который существует в цикле Д алгоритма Карлсона на стр. 286.28613.9 Редуцированный фильтр Бирмана13.9Редуцированный фильтр БирманаВ приложениях выбор измерительных средств бывает ограничен так, чтов измерение z попадает лишь часть (заранее известная) элементов оцениваемого вектора x.
Пусть эта часть — первые qj < n элементов для j-й строкиматрицы наблюдений:hTj = [ |? ?{z· · · ?} 0 · · · 0 ] = [ ? ? · · · ? 0 · · · 0 ].(13.9)qj <n , hqj 6=0Теорема 13.6 ([78]).Пусть калмановская процедура скалярногоbDbLbT и Pe = LeDeLeT ,обновления (13.5), (13.6) использует разложения Pb = Lгде L — нижние треугольные с единичной диагональю матрицы, D — диагональные (с положительными элементами) матрицы. Тогда данная процедура(13.5), (13.6) эквивалентна пунктам В, Г, Д и Е, следующего алгоритма.II. Обработка измерения:e = D(t− ), xe = L(t− ), De = x(t−А.
Начальное присваивание: Lii ).iБ. m-кратное повторение процедуры скалярного обновления.Для j = 1, 2, . . . , m выполнять:eT h ; v = Dfe . Их вид:В. Вычислить векторыf =Lf = [ ? ? · · · ? 0 · · · 0 ]T , v = [ ? ? · · · ? 0 · · · 0 ]T .Г. Задать начальные значения α 0 = r; K = [0 . . . 0]T .Д. Для i = q + 1, q + 2, . . .
, n выполнять: Ki := eliq vq .Е. Для i = q, q − 1, . . . , 2, 1 выполнять:началоα := α 0 + vifi ; dbi := dei α 0 /α ;λ := −fi/α ;Ki := vi ;bli := eli + λK ; K := K + eli v i ;(✮)0α := α .конецЖ. Вычислить векторы ν := (z − hT xe)/α ; xb := xe + Kνс экстраполяцией между повторениями:e := Lb; De := Db; xLe := xb.З. Завершающее присваивание:bbL(t+D(t+xb(t+b.i ) := L ;i ) := D ;i ) := x28713 Устойчивые алгоритмы фильтрацииЗдесь hT — j-я строка матрицы H(ti ), имеющая вид (13.9); z —j-й элемент вектора z(ti ); r — j-й элемент rj (ti ) диагональной матрицы ковариаций шума измерений R(ti), j = 1, 2, . . .
, m — номерскалярного измерения в составе вектора измерений z(ti ) в моментвремени ti .Замечание 13.7. Строка ✮ в алгоритме на стр. 287 выглядит также, как в замечании 13.5 на стр. 282.Замечание 13.8. Представление (13.9) мотивирует использованиеLD-версии фильтра Бирмана по подразд. 13.7. Если жеhTj = [ 0 · · · 0 ? ? · · · ? ],то сокращенный объем вычислений будет получаться при использованииU D-версии фильтра Бирмана [15].Замечание 13.9. Сокращение объема вычислений в редуцированных таким образом версиях получается наибольшим среди алгоритмов этогокласса, поскольку значение q берется как qi — индивидуально для каждогоиз m повторений цикла Е в алгоритме на стр.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
