Вычислительные методы алгебры и оценивания. И.В. Семушкин (2011) (1185350), страница 49
Текст из файла (страница 49)
. . , fn ) ,.T.K1 = Ũ11f1 . 0 . . . 0 .Для j = 1, . . . , n выполнятьαj = αj−1 + fj2 ,βj = (αj−1 /αj )1/2,γj = fj /(βj αj ) ,Ûjj = βj Ũjj ,Ûij = βj Ũij − γj Kj (i),i = 1, 2, . . . , j − 1, j 6= 1 ,Kj+1(i) = Kj (i) + fj Ũij ,i = 1, 2, . . . , j .По завершении цикла выполнитьK = Kn+1/αn ,x̂ := x̃ + K(z − aT x̃) .31113 Устойчивые алгоритмы фильтрацииВариант 7. Факторизованный LLT ковариационный алгоритмКарлсона. Найдите его на стр. 285, подразд.
13.8.Вариант 8. Редуцированный фильтр Бирмана. Найдите его на стр. 287,подразд. 13.9.Вариант 9. Редуцированный стандартный ковариационный алгоритмБар-Ицхака–Калмана. Найдите его на стр. 288, подразд. 13.10.Замечание 13.23. В этом варианте используются только способы 2и 4 генерации матрицы A (см. выше п. 2 в подразд.
13.14). Инициализацияи обработка первых n наблюдений выполняются, как в варианте 1. Послеэтого hкаждое очередноенаблюдениеимеет iвид z = (a(q) )T x̄(q) + v, так какihTTTaT = a(q) ... 0 и x̄T = x̄(q) ... x̄(s) , где s = n − q. Вектор x̄, также как и соответствующие ему векторы x, x̃, x̂, разбит на две части вида:x̄(q) размерности q и x̄(s) размерности s. Запишем алгоритм Бар-Ицхака–Калмана для варианта 9:T(ii) Обработка наблюдений (очередные данные) z = a(q) x̄(q) + v:Tα = a(q) P̃ (qq) a(q) + r,K (q) = P̃ qq a(q) /α ,TP̂ (qq) = P̃ (qq) − K (q) a(q) P̃ (q) ,Tx̂(q) = x̃(q) + K (q) (z − a(q) x̃(q) ) ,K (sq) = P̃ (sq) (P̃ (qq) )−1 ,P̂ (sq) = K (sq) P̂ (qq) ,P̂ (ss) = P̃ (ss) − K (sq)(P̃ (qq) − P̂ (qq) )(K (sq))T ,x̂(s) = x̃(s) + K (sq) (x̂(q) − x̃(q) ) .Здесь одинарный верхний индекс в скобках указывает размерность вектора,а двойной верхний индекс в скобках указывает размер матрицы.
МатрицыP̂ и P̃ из варианта 1 здесь разбиты на блоки по схеме: (qq) (qs) PPP =,P = PT.(sq)(ss)PPВариант 10. Редуцированный стабилизированный ковариационныйалгоритм (Бар-Ицхака–Джозефа). Как и в варианте 9, здесь применяетсядекомпозиция, т. е. разбиение на блоки векторов и матриц, что выделяет31213.15 Варианты задания на лабораторный проект № 9редуцированную часть алгоритма вида (13.12) и специфическую часть, наподобие (13.13).
Для инициализации и обработки первых n измерений используется алгоритм варианта 2.Вариант 11. Редуцированный квадратно-корневой алгоритм БарИцхака–Поттера, с нижнетреугольным разложением, S ≡ L. За основуберется алгоритм варианта 9, но в нем редуцированная стандартная частьзаменяется на редуцированную часть размера (q × q). Формулы специфической части заменяются, соответственно, на следующие выражения:K (sq) = L̃(sq) (L̃(qq) )−1 ,L̂(sq) = K (sq) L̂(qq) ,L̂(ss) = L̃(ss) .Этап экстраполяции выполняется, как в алгоритме Поттера, вариант 3, приэтом (qq)L0S≡L=L(sq) L(ss)и L(qq) , L(ss) — нижние треугольные матрицы (с верхним знаком ˜ или ˆ).Для инициализации и обработки первых n измерений используется алгоритмварианта 3.Вариант 12.
Редуцированный квадратно-корневой алгоритм (БарИцхака–Бирмана–Медана), с P = LDLT -разложением по типу алгоритмаБирмана (см. вариант 4).Здесь используются декомпозиции вида: (qq) (qq) qP(P (sq) )TL0D 0P =, L=, D=,0 DsP (sq) P (ss)L(sq) L(ss)где L — нижние треугольные с единичной диагональю матрицы, D — диагональные матрицы. Действует замечание 13.23 после Варианта 9, но дляинициализации и обработки первых n наблюдений используется алгоритмварианта 4. Для обработки каждого из остальных измерений применяетсяследующий алгоритм:(1) Выполнить алгоритм из варианта 4, но в применении к векторам размерности q и матрицам размера (q × q), т.
е. по данным z, a(q) , L̃(qq) ,D̃(q) , x̃(q) найти L̂(qq) , D̂(q) , x̂(q) .31313 Устойчивые алгоритмы фильтрации(2) ВычислитьK (sq) = L̃(sq) (L̃(qq) )−1 ,L̂(sq) = K (sq) L̂(qq) ,L̂(ss) = L̃(ss) ,D̂(s) = D̃(s) ,x̂(s) = x̃(s) + K (sq) (x̂(q) − x̃(q) ) .Вариант 13. Редуцированный квадратно-корневой ковариационный алгоритм (Бар-Ицхака–Карлсона). Действует замечание 13.23 после Варианта 9, но для инициализации и обработки первых n измерений используетсяалгоритм варианта 7.
Для обработки остальных измерений применяется следующий алгоритм:(1) Выполнить алгоритм из варианта 7 (его следует вывести самостоятельно, наподобие алгоритма варианта 6), но в применении к векторамразмерности q и матрицам размера (q × q), т. е. по данным z, a(q) , L̃(qq) ,x̃(q) найти L̂(qq) и x̂(q) .(2) Вычислить все остальные матрицы и вектор x̂(s) .Вариант 14. Стандартный информационный алгоритм (см. стр. 245).(i) Инициализация (начальные значения x0, Λ0 ) :....d0 = Λ0 x0 ,Λ . d := Λ0 . d0 .(ii) Обработка наблюдений (очередные данные z = aT x̄ + v):....TΛ . d := Λ .. d + a a ..
z /r .(iii) Выдача результата: x̂ = Λ−1 d .В качестве начальных значений рекомендуется взять x0 = 0, Λ0 = 0.Вариант 15. Квадратно-корневой информационный алгоритм. См. подразд. 7.3, стр. 110, формулу (7.7), здесь — ее рекуррентная версия (13.28).(i) Инициализация (начальные значения R̃0, x0) :hi hiz̃0 = R̃0 x0;R̂0ẑ0 = R̃0z̃0 .31413.15 Варианты задания на лабораторный проект № 9(ii) Обработка наблюдений (очередные скалярные данные z = aT x̄ + v):R̂j ẑjR̂j−1 ẑj−1= Tj, j = 1, 2, . .
. , m ,(13.28)0 eaTzгде Tj — ортогональная матрица, которая выбирается так, при чтобыR̂j−1,j—вести к верхнетреугольному виду расширенную матрицуaTномер очередного измерения, все матрицы R здесь имеют размер (n×n),при этом все Rj , j ≥ 1, — верхние треугольные.(iii) Выдача результата: x̂ = R̂j−1 ẑj .Hачальные значения рекомендуется взять в виде x0 = 0, R̃0 = 0. В качестве преобразований Tj возможны четыре процедуры, указанные в описанииварианта 3.
Таким образом, всего имеем 4 разновидности данного варианта:TjвариантХаусхолдер15.1Гивенс15.2ГШО15.3МГШО15.4Замечание 13.24. В каждом варианте лабораторного проекта № 9необходимо выполнить проверку всех положений для задачи с мультиколлинеарностью, которые изложены в подразд. 13.13 на стр. 300–305.31514Ортогонализованные блочныеалгоритмыВ этом разделе представлены дальнейшие инновации в области вычислительных методов оценивания, основанные на соединении идей факторизации матриц, скаляризации процесса обработки данных и блочной ортогонализации массивов данных.
В оригинале эти алгоритмы называются arrayalgorithms [17]. Использован материал диссертации М. В. Куликовой [51] ссистематизацией (переименованиями) некоторых алгоритмов.14.1Задача оцениванияПерепишем постановку задачи оценивания, определенную общими выражениями (13.1), с указанием дискретного времени не в скобках, а в индексе.Пусть дискретная динамическая система описывается уравнениями:xt+1 = Φt xt + Bt ut + Gt wt ,(14.1)zt = Ht xt + vt,(14.2)где xt — n-мерный вектор состояния системы, zt — доступный m-мерныйвектор измерений, ut — детерминистский r-мерный вектор (входное управляющее воздействие), t = 0, 1, .
. . , — дискретное время. Матрицы Φt ∈ Rn×n ,Bt ∈ Rn×r , Gt ∈ Rn×q , Ht ∈ Rm×n , Qt ∈ Rq×q , Qt ≥ 0 и Rt ∈ Rm×m ,Rt > 0 известны, могут быть параметризованы и полностью характеризуют систему (14.1), (14.2). Последовательности {w0, w1, . . . } и {v1, v2, . . . } —независимые нормально распределенные последовательности шумов с нулевыми средними значениями и ковариационными матрицами Qt и Rt , соответственно.
Без потери общности считаем, что {w0, w1, . . . } и {v1, v2, . . . }не зависят от начального вектора состояния системы x0, распределенного14.2 Блочные алгоритмы в исторической перспективепо нормальному закону с математическим ожиданием x̄0 и ковариацией P0 ,т. е. x0 ∼ N (x̄0, P0 ).Задача оценивания заключается в получении оценки неизвестного вектора xt, t = 0, 1, . . . по доступным наблюдениям Z1N = {zt , t = 1, 2, . . . , N },содержащим информацию о векторе xt . Если N = t − 1, оценка x̂t вектора называется экстраполяционной оценкой, точнее, оценкой одношаговогопредсказания. Если N = t, оценка x̂t вектора называется отфильтрованнойоценкой.
Если N > t, оценка x̂t вектора называется сглаженной оценкой.Задачи сглаживания в этом пособии не рассматриваются. Они подробнопредставлены в [61], где можно видеть, что в составе алгоритмов сглаживания присутствуют оптимальные оценки от фильтра Калмана.Оптимальные оценки — это те, которые наилучшим образом (в заранее определенном смысле) соответствуют истинному значению вектора xt .Если критерий качества оценивателяусловным математиче T N Jt определенским ожиданием как Jt , E x̃t x̃t Z1 = E (xt − x̂t )T (xt − x̂t ) Z1N , гдеE {·} — оператор математического ожидания, x̃t , xt − x̂t — погрешность, тооцениватель x̂t , минимизирующий критерий Jt , называется оптимальным всреднеквадратическом смысле оценивателем.14.2Блочные алгоритмы в исторической перспективеДискретный стандартный ковариационный фильтр (СКФ) Калманадается уравнениями (13.3).
Их вывод широко известен, см. например,[61, 113]. В обозначениях системы (14.1), (14.2) эти уравнения выглядятнесколько более обозримо в следующей записи:•Этап экстраполяции: t = 0, 1, . . . ; P0+ = P0 ,+оценка: x̂−t+1 = Φt x̂t + Bt ut ,−ковариация ошибки: Pt+1= Φt Pt+ ΦTt + Gt QtGTt .•(14.3)(14.4)Этап обработки измерений (фильтрация): t = 1, 2, . . . ,Kt = Pt− HtT (HtPt− HtT + Rt )−1,−−оценка: x̂+t = x̂t + Kt (zt − Ht x̂t ),ковариация: Pt+ = Pt− − Kt Ht Pt− ,(14.5)(14.6)(14.7)где все данные берутся из модели (14.1), (14.2): Φt , Gt , Bt , Ht , Qt ,Rt — известные матрицы-параметры линейной дискретной динамической31714 Ортогонализованные блочные алгоритмысистемы, ut — вектор управления, zt — доступный вектор наблюдений, xt —вектор состояния системы с начальным значением x0 ∼ N (x̄0, P0 ).Соединяя в (14.3)–(14.7) два этапа в один, для t = 0, 1, .