О современных ортогонализованных алгоритмах оптимальной дискретной фильтрации Текст научной статьи по специальности «Математика»

УДК 519.61+519.711:681.5 Б01: 10.14529/ттр180401

О СОВРЕМЕННЫХ ОРТОГОНАЛИЗОВАННЫХ АЛГОРИТМАХ ОПТИМАЛЬНОЙ ДИСКРЕТНОЙ ФИЛЬТРАЦИИ

Ю.В. Цыганова1, М.В. Куликова2

Ульяновский государственный университет, г. Ульяновск, Российская Федерация 2

В настоящее время вычислительные методы оптимального оценивания стали самостоятельной областью исследования и получили большое развитие. Современные численно эффективные ортогонализованные алгоритмы привлекательны не только своей устойчивостью к ошибкам машинного округления, но и приспособленностью алгоритмов, использующих различные типы матричных ортогональных преобразований, к программной реализации на параллельных или векторных вычислительных системах. Эти свойства позволяют разрабатывать новые эффективные информационные технологии, в частности, при решении задач в режиме реального времени и при обработке больших данных. Статья содержит краткий обзор современных ортого-нализованных алгоритмов оптимальной линейной дискретной фильтрации. Рассмотрены четыре класса ортогонализованных алгоритмов: квадратно-корневые ортогонализованные алгоритмы, алгоритмы на основе методов взвешенной ортогонализации, 7-ортогонализованные алгоритмы и алгоритмы на базе методов сингулярного разложения. Приведена классификация алгоритмов по типам матричных ортогональных преобразований, на основе которых эти алгоритмы построены. Такая классификация позволяет легче понять метод построения ортогонализованного фильтра и выбрать способ эффективной программной реализации при решении практических задач в классе многомерных дискретных линейных стохастических систем. В работе исследованы вычислительные аспекты ортогонализованных алгоритмов: численная устойчивость к ошибкам машинного округления и способы эффективной программной реализации. Все рассмотренные алгоритмы являются алгебраически эквивалентными стандартной реализации дискретного фильтра Калмана, но существенно превосходят его по своим вычислительным свойствам. Полученные результаты сравнительного исследования позволяют сделать вывод о том, что применение ортогонализованных алгоритмов при решении практических задач помогает получить численно эффективные решения.

Ключевые слова: дискретная фильтрация; линейные стохастические системы; фильтр Калмана; матричные ортогональные преобразования; ортогонализованные алгоритмы.

Введение

В любой области научных исследований, связанной с математическим моделиро-

ВсШИбМ^ Д&ННЫб^

собираемые в результате наблюдения за некоторым объектом, служат исходным материалом для решения таких задач, как построение математической модели, основанное на модели оценивание состояния и основанное на оценках управление объектом.

Поскольку

данные всегда первичны^ для любой науки, имеющей дело с данными, моделирование является одной из главных задач. Модель извлекает информацию из данных благодаря тому, что она выявляет отношения между измеряемыми

характеристиками объекта и тем самым позволяет (если она не дескриптивная, а пре-диктивная), предсказывать свойства или поведение объекта всякий раз до того, как о нем поступят новые результаты измерений. Предиктивная модель, как правило, выражает информацию об объекте в математической форме, поэтому использование таких моделей справедливо относят к области математических информационных технологий. Данная работа посвящена современным методам оптимальной дискретной фильтрации в задаче оценивания состояний многомерных дискретных линейных стохастических моделей сложных систем на основе измерительных данных.

Фильтр Калмана (ФК) [1] является наиболее известным математическим рекуррентным алгоритмом оценивания вектора состояния линейной стохастической системы на основе доступных измерений. В настоящее время существует широкий класс робастных (численно устойчивых к ошибкам машинного округления) модификаций ФК, которые основаны на различных методах факторизации ковариационных матриц. Большое количество научных публикаций в области параметрического оценивания подтверждает практическую эффективность данных методов.

Среди всех робастных модификаций ФК можно выделить класс ортогонализован-ных алгоритмов, которые превосходят по своим вычислительным свойствам стандартные алгоритмы калмановской фильтрации. Современные ортогонализованные алгоритмы обладают улучшенными вычислительными свойствами, включая устойчивость к ошибкам машинного округления, компактную и однородную структуру, которая непосредственно ориентирована на единообразную обработку всех д&нных и идеальна для программной реализации средствами параллельного программировал ния^ что позволяет эффективно обрабатывать большие объемы ДсШ н ых.

Данная ст&тья, являясь ^ по существу^ обзорной, имеет в своей основе изложение тех ключевых математических идей, которые привели к самым эффективным вычислительным алгоритмам во всех приложениях, где так или иначе возникает потребность в применении методов оптимальной дискретной фильтрации.

Статья организована следующим образом. В п. 1 излагается задача оптимальной дискретной фильтрации состояния многомерной дискретной линейной стохастической системы и рассматриваются основные типы алгоритмов дискретной фильтрации. В п. 2 сформулированы основные идеи, позволяющие преодолеть численную неустойчивость стандартного фильтра Калмана. В п. 3 обсуждаются вопросы устойчивости и эффективной программной реализации ортогонализованных алгоритмов.

1. Предпосылки

Рассмотрим многомерную дискретную линейную стохастическую систему 6, представленную в пространстве состояний уравнениями

где к - дискретный отсчет времени, хк Е Кга - вектор состояния, подлежащий оцениванию; хк Е - доступный вектор измерений; ик Е - детерминированное управление. Шумы } и [ук} - независимые нормально распределенные последовательности с нулевым математическим ожиданием и матрицами ковариаций Qk > 0 и Як > 0, соответственно. Начальное значение вектора состояния Хо является гауссовым слу-

Xk+i = Fk xk + Вк Uk + Gk Wk, zk+1 = Hk+ixk+1 + vk+i, к > 0,

(1) (2)

чайным вектором с характеристиками х0 ~ N(х0, По). Вектор х0 те зависит от {ши} и {ьи}■

Дискретная линейная стохастическая система 6 (1), (2) является математической моделью объектов различной природы. Поскольку значения вектора состояния недоступны прямому наблюдению ввиду присутствия помех в (2), мы можем вычислить только оценки вектора состояния в результате решения задачи оценивания (или фильтрации). Для линейных систем минимум квадратичного критерия качества оценивания соответствует оптимальной в среднеквадратическом смысле оценке Хи вектора состояния хи, которую можно вычислить с помощью известного с середины прошлого века алгоритма Р. Калмана.

1.1. Стандартные алгоритмы дискретного фильтра Калмана

Свое знаменитое решение задачи оптимальной дискретной фильтрации Рудольф Эмиль Калман [1] предложил в начале 1960-х годов прошлого века. Это решение получило название <фильтр Калмана >. Изначально фильтр Калмана был построен для решения задачи параметрического оценивания движущихся технических объектов. И хотя его первое практическое использование было СВЯЗсШО с областью аэронавтики, впоследствии данный метод нашел свое применение в самых разных областях науки и техники, таких как обработка изображений, спутниковая геодезия, телекоммуникационные сети, метеорология, эконометрика, молекулярная биология и во многих других. В наше время калмановская фильтрация является одним из самых популярных методов решения различных прикладных задач, связанных с параметрическим оцениванием по экспериментальным данным.

Стандартный фильтр Калмана записан в так называемой ковариационной форме, поскольку на каждом шаге в алгоритме рекуррентно обновляются значения коварна-ционных м&триц ошибок оценивания вместе с обновлением оценок вектора состояния. Алгоритм состоит из двух этапов: экстраполяция1 и фильтрация2. Эту форму ФК называют двухстадийной. Соответствующие уравнения приведены ниже.

CKF-TM - Conventional Time and Measurement Update Kalman Filter

® Начальные данные

положить Жо|о = xo, Po|o = По (По > 0)

©Экстраполяция (k > 0)

оценка ковариация ошибки Xfc+i = Fk Xk\k + Bk uk Pk+i = Fk pkik Fk + Gk Qk G~[

©Фильтрация (k > 0)

коэффициент Калмана Kk+1 = Pk+lH]^+1Re,1+v Re,k+1 = Rk+1 + Hk+1Pk+1Hk+1

оценка ковариация ошибки Xk+i|k+i = Xk+1 + Kk+1ek+b ek+1 = zk+1 — Hk+1Xk+1 pk+1|k+1 = pk+1 — Kk+1Hk+1pk+1

Условия сходимости алгоритма Калмана подробно рассмотрены в [14].

Алгоритм СКК-ТМ можно представить в сокращенной одностадийной форме [20 р. 895], которая требует меньше вычислений, чем двухстадийная форма алгоритма.

1 Предсказание значений оценки вектора состояния и матрицы ковариации ошибки оценки по предыдущим доступным измерениям.

2 Обновление предсказанной оценки вектора состояния и матрицы ковариации ошибки оценки в соответствии с полученным в текущий момент времени измерением.

Алгоритмы фильтрации информационного типа вместо матрицы ковариации ошибки оценивания Pk на каждом шаге рекуррентно обновляют значения информационной матрицы где Лк = P-1. Вместо оценки вектора состояния вычисляют информационную оценку dk = Л kxk. В отличие от алгоритмов ковариационного типа, информационные алгоритмы могут стартовать в отсутствии априорной информации, т.е. при Л0 = 0 [23, р. 262].

Быстрые алгоритмы играют важную роль в задаче оценивания стационарных линейных стохастических систем большой размерности, т.е. при n >> 0. В их основе лежит следующий подход. Поскольку для стационарных систем ранг матрицы 8Pk = Pk+i — Pk часто много меньше, чем ранг Pk, то на каждой итерации алгоритма целесообразно обновлять значения матрицы 6Pk вместо Pk.

Согласно [22, лемма 11.1.1], для эрмитовой матрицы 8P0 ранга а существует разложение вида SP0 = P1 — П0 = L0M0 Lq, где M0 Е Raxa - эрмитова невырожденная матрица. Тогда матрицы 6Pk можно представить в виде 6Pk = LkMkLj, к > 0, где Mk Е Raxa - эрмитова невырожденная матрица. Для Lk и Mk справедливы следующие рекуррентные соотношения (к > 0): Lk+1 = (F — KkR-^H)Lk, Mfc+1 = Mk — Mk LjHTR-l+1 HLk Mk.

Отметим, что в классе быстрых алгоритмов разложение SP0 = L0M0LJ осутдвств-

L0 M0

разностного уравнения Риккати уже не требуется, поскольку в алгоритме обновляются матрицы меньшего размера Lk, Mk, где Lk Е Mrax^ Mk Е Maxa. По ним всегда можно восстановить матрицу ковариации Pk+1 для любо го к = 1,... ,K. В общем случае быстрый алгоритм фильтрации задается рекуррентными уравнениями [22, гл. 11], где показано, что в случае m < n, q < n, который часто встречается на практике, сложность алгоритма фильтрации сокращается с O(n3^о O(an2). Как правило, а << n для систем большой размерности.

1.2. Численная неустойчивость фильтра Калмана

Фильтр Калмана практически сразу же после своего появления завоевал большую популярность среди специалистов в области параметрического оценивания, поскольку представлял собой вычислительный алгоритм, идеальный для реализации на ЭВМ. Однако, он имеет один существенный недостаток, а именно, неустойчивость к ошибкам машинного округления. В теории вычислительные свойства фильтра Калмана характеризуются ковариационной матрицей ошибки оценивания, которая является решением матричного дифференциального или разностного уравнения Риккати. С момента первых реализаций фильтра Калмана на ЭВМ сразу же было обнаружено, что вычисленные среднеквадратические ошибки оценивания часто существенно отличаются от теоретических значений. В некоторых случаях полученное решение уравнения Риккати имело отрицательные диагональные элементы, что теоретически невозможно. Проблема заключалась в ошибках машинного округления.

Данная проблема хорошо известна специалистам в области оптимальной дискретной фильтрации. Теоретический численный анализ фильтра Калмана проведен в из-В6СТНОИ работе [15]. Во многих работах на практических примерах показана численная неустойчивость ФК, см. например, подробный анализ в [23, р. 205, example 6.2]. Обсуждение указанной проблемы можно найти в [6,13,17,22,23,28].

3Данное разложение можно найти с помощью алгоритма Bunch-Kaufman [24].

Ошибки машинного округления, которые неизбежно возникают из-за ограниченности разрядной сетки при реализации арифметики с плавающей запятой, характеризуются единичной ошибкой округления е0кр4- При реализации алгоритма Калмана на ЭВМ был обнаружены проблемы вычислении ^ связанные с ошибками машинного округления. Для того, чтобы программная реализация ФК давала правильный результат, теоретическое решение уравнения Риккати, согласно которому обновляется матрица ковариаций ошибок оценивания Ри, должно совпадать с вычисленным на ЭВМ значением. Если это условие не выполняется, соответствующая задача является плохо обусловленной. В [23] авторы приводят несколько причин численной неустойчивости ФК, среди которых неопределенности в параметрах системы (ошибки моделирования), большой диапазон фактических значений параметров системы вследствие неудачно построенной математической модели, плохая обусловленность ковариаци-ониои матрицы невязки измерении ^ низкая точность представления вещественных чисел в ЭВМ и др. В [17] указан ряд причин потери точности при программной реализации ФК, в том числе потеря свойства положительной определенности ковариационными матрицами ошибок оценивания, плохая обусловленность матриц, требующих операции обращения, при точных измерениях, а также при проведении в относительно короткое время ряда неточных измерений с примерно одинаковым вкладом их в пересчет оценки вектора состояния.

Перечисленные выше причины нередко появляются при решении практических задач. Различия между теоретическими и вычисленными значениями в алгоритме Калмана зачастую становятся настолько существенными, что возникает вопрос о целесообразности использования полученных результатов.

К сожалению, окончательно избавиться от проблемы ошибок машинного округления в настоящее время еще невозможно ввиду ограниченности разрядной сетки компьютера. Но можно существенно уменьшить их влияние с помощью численно эффективных модификаций фильтра Калмана. Численная эффективность означает, в первую очередь, устойчивость по отношению к ошибкам машинного округления, сокращение объема вычислений, а также наличие эффективных (компактных и одно-р од н ых) структур данных, поступающих для обработки в алгоритм фильтрации.

Первые численно эффективные модификации ФК появились уже вскоре после его открытия. Основной идеей построения таких алгоритмов является применение методов факторизации ковариационных матриц ошибок оценок, участвующих в уравнениях фильтра. Первым так и м ал го ритм о м был квадратно-корневой алгоритм Пот-тера [2], который принято считать отправной точкой для развития целого направления в области линейного оценивания и управления, которое продолжает активно развиваться и сейчас. Данный факт подтверждается большим количеством научных публикаций. Отметим хорошо известные авторам монографии, содержащие описание численно эффективных модификаций ФК: зарубежные [13, 22, 23, 28] и отечественные [17,29,30]. Основные направления в теории устойчивых методов оптимальной дискретной фильтрации приведены в [23] и на русском языке в недавней работе [31].

Численно эффективные реализации фильтра Калмана широко используются в мире для решения практических задач, о чем свидетельствует большое количество

4Величина еокр определяется как наибольшее число такое, что с учетом машинного округления в операциях с вещественными числами 1 + еокр = 1 или 1 + £0кР/2 = 1. Значение этой ошибки зависит от конкретного способа реализации машинной арифметики.

зарубежных научных публикаций, в том числе за последние годы, например, управление движением мобильных роботов [26], задачи слежения за техническими объектами [27], обработка цифровых сигналов [35], оценивание волатильности курсов валют [32], задачи судовождения [33]. В то же время в данной области исследовании отечественных публикаций значительно меньше. Авторы надеются ^ что данная статья сможет в некоторой степени восполнить образовавшийся пробел о мировых достижениях в области численно эффективных методов оптимальной дискретной фильтрации.

2. Современные ортогонализованные алгоритмы

В настоящее время наиболее предпочтительными для решения практических задач являются численно эффективные ортогонализованные алгоритмы [22].

Согласно [37], ортогонализованным алгоритмом назовем «алгебраически эквивалентный дискретному фильтру Калмана алгоритм вычисления оценок хк и/или хк\к вектора состояния дискретной линейной стохастической системы (1), (2) и ковариационной матрицы ошибки оценивания Рк и/ил и Рк\к, основанный на матрично-блочном представлении разностного уравнения Риккати и матричном ортогональном преобразовании как способе вычисления указанных величин >.

Как и все известные устойчивые модификации дискретного фильтра Калмана, ортогонализованные алгоритмы фильтрации обладают двойной точностью вычислений, а использование численно устойчивых ортогональных преобразований на каждом шаге дает дополнительные преимущества6. Кроме того, однородная структура современных ортогонализованных алгоритмов идеально подходит для распараллеливания на многопроцессорных ЭВМ [16,22,35].

В ортогонализованных алгоритмах вместо ковариационной матрицы ошибки оценивания Рк вычисляют матричные факторы разного вида. Например, квадратно-корневые ортогонализованные алгоритмы основаны на разложении Холецкого матрицы ковариации Рк = Брк БРк и на матричных ортогональных преобразованиях вида QA = Я (либо QA = Ь) как способе обновления матричных факторов БРк7. Такая организация вычислений сохраняет положительные элементы на диагонали ковари-Рк

обеспечивают вдвое больше точных цифр после десятичной запятой, чем стандартная реализация фильтра Калмана. Увеличение точности вычислений происходит за счет того, что в ортогонализованных алгоритмах применяются ортогональные преобразования для обновления матричных факторов БРк в соответствии с уравнением Риккати, записанным в блочно-матричной форме

Q

Фактор Холецкого БРк вычисляют с помощью любого (^-преобразования блочной прямоугольной матрицы. Такая форма матричных уравнений с учетом различных размерностей позволяет строить разные варианты алгоритмов.

" A ' SPk

B 0

5В иностранной литературе используют термин «array algorithms for Kalman filtering».

6

Подробное обсуждение этого вопроса можно найти в [22].

7Разложение Холецкого выполняется только один раз для начальных значений Р0 > 0, и далее

о

обновление уравнений фильтра происходит в терминах матричных квадратных корней Яр

к '

Первым ортогонализованным алгоритмом принято считать квадратно-корневой информационный фильтр SEIF (Square Root Information Filter) [4]. Ортогонализован-ные алгоритмы ковариационного и смешанного типов были предложены и развиты Т. Кайлатом и соавторами [20].

Организация вычислении в блочной ортогонализованной форме дает алгоритму следующие преимущества: 1) позволяет обрабатывать данные с двойной точностью; 2) на каждом шаге рекурсии используются численно устойчивые ортогональные преобразования, что гарантирует сохранение специального верхнего/нижнего треугольного вида матриц SPk = рЦ2'-, 3) данные методы свободны от операции матричного обращения, чувствительной к ошибкам округления.

Далее рассмотрим современные классы ортогонализованных алгоритмов.

2.1. Квадратно-корневые ортогонализованные алгоритмы

В начале прошлого века А.Л. Холецкий [23, р. 373] предложил оригинальное решение задачи МНК с помощью метода извлечения квадратных корней из неотрицательно-определенной симметрической матрицы, позволяющего представить ее в виде произведения нижней треугольной матрицы Ь и верхней треугольной матрицы и = Ьт. В свою очередь, методы матричных ортогональных преобразований первоначально были разработаны также для решения задачи МНК. В калмановской фильтрации данные методы применяются для вычисления факторов Холецкого ковариационных матриц ошибок оценивания на этапах экстраполяции и фильтрации квадратно-корневого алгоритма.

С. Шмидт [5] впервые предложил использовать метод Холецкого для вычисления матричного квадратного корня ¿^такого, что Брк Брк = Рк. Есл и Брк1к - фактор Холецкого матрицы Рк\к, а БQk - фактор Холецкого матрицы ковариации шума Qk, тогда на этапе экстраполяции уравнение для матрицы Рк+\ можно переписать в виде Брk+1Бр'k+1 = ^к5рк\к5р'к\к^к + ^иБQk5~о;к■ Метод вычисления матрицы Брк, пред-ложенныи Шмидтом, заключается в ортогональном преобразовании прямоугольной (2 х 1)-блочной матрицы вида

Q

S~Pk\k Fk

qT riT

SQk Gk

SPTk+i 0

(3)

где Q - матрица ортогонального преобразования. Шмидт предложил алгоритм для построения матрицы Q, получивший название модифицированный алгоритм Грама - Шмидта [23].

П. Каминский [7] предложил новые модификации ковариационных и информационных квадратно-корневых алгоритмов с улучшенными вычислительными свойствами. Учитывая свойство двойственности ковариационных и информационных алгоритмов, он разработал новую вычислительную схему

Q

SRk 0

H к Spk Spk

SRe,k

Kf,k

0

Spk

(4)

для этана обработки измерений в дискретном фильтре, где Бре к — квадратный корень

матрицы Яе,к ковариации невязки измерений ek, Kf,k = РкН^ Яе1/2 и Q матрица ортогонального преобразования к нижнему треугольному виду блочной матрицы в

левой части (4). Он также предложил использовать на этапе экстраполяции выражение (3).

Известно, что основной объем вычислении в фильтре Калмана приходится на решение разностного уравнения Риккати. Для его сокращения М. Морф и Т. Кай-лат [10] выдвинули идею совмещения в одном ортогонализованном алгоритме этапов экстраполяции и фильтрации. Они предложили вычислять фактор Холецкого БРк+1 по известному фактору Холецкого БРк, полученному на предыдущем шаге алгоритма, с помощью процедуры триангуляризации блочной матрицы

Г oT 0 ST e,k

Q T у S Tk HT STk Fk = 0

0 cT лТ SQk Gk . 0

Rlk

SPk+i

(5)

где Q - матрица ортогонального преобразования, приводящего к верхнему треугольному виду матрицу в левой части (5), при этом КС^ к = ГкКfкк.

Изложенные идеи касаются алгоритмов ковариационного типа. Идея квадратно-корневого алгоритма информационного типа принадлежит П. Дайеру и С. МакРей-нольдсу [4]. Они построили алгоритмы экстраполяции и фильтрации для информационного фильтра в терминах матричных квадратных корней. Алгоритм ЭШР, получивший широкое распространение также, как и квадратно-корневые ковариационные фильтры, обладает улучшенными вычислительными свойствами по сравнению со стандартной реализацией информационного фильтра. Уравнения совмещенного (одностадийного) квадратно-корневого ковариационного алгоритма представлены ниже.

SRCF - Square-Root Covariance Filter

положить вычислить

Х0 = Й0, Po = По > 0 (R > 0 Qo > 0)

Sn0, Sr0, Sq0 - нижние треугольные факторы Холецкого

^^ Рекуррентное обновленрте данных (к > 0)

Г sTk 0 ] ST e,k

выполнить QkAk = Rk Qk STk Hk FkT 0 STk gT - — 0 0 s Pk+1 0

оценка xk+i = Fk xk + Bk Uk + Kp,k ek, ek = S-ik (zk - Hkxk)

вычислить STk+u SQk+i

Алгоритм SRCF можно сформулировать в двухстадийной ковариационной форме в виде алгоритма SECF-TM (Time and Measurement Update SRCF) [22, глава 12]. Сравнение вычислительной сложности квадратно-корневых алгоритмов см. в [15].

0

2.1.1. Расширенные ортогонализованные алгоритмы

Расширенные алгоритмы являются наиболее подходящими для решения задач на современных вычислительных комплексах с параллельной архитектурой [16,22]. Основная идея таких реализаций ФК заключается в формировании пред-массива (рге-аггау), который представляет собой блочную матрицу, и последующем его приведении к требуемому нижнему или верхнему (либо их комбинации) треугольному блочному пост-массиву (post-array) с помощью матричных ортогональных преобразований. В работе [20,22] Т. Кайлат и соавторы предложили несколько вариантов расширенных

ортогонализованных алгоритмов. Уравнения расширенного квадратно-корневого одностадийного ковариационного фильтра (алгоритм еЭКСР) представлены ниже.

eSRCF - extended Square-Root Covariance Filter

положить вычислить найти нач. оценку

Хо = xo, Po = По > 0 Ro > 0 Qo > о

Sn0, Sr0, Sq0 - нижние треугольные факторы Холецкого

¿0 = sn0lx0

^^ рекуррентное öbhöbjiehpie данных k > 0^

s—kZk SRe,k 4k -ek+1

Zk = 0 с t SPk+i Zk+1

0 0 0 Yk

выполнить qkak = rk

Qk

Sjk.

0

S~pk Hk sp Fkp 0 spkGpp

вычислить

SRk+n SQk+!

В алгоритме eSRCF в каждый момент времени к доступна оценка вектора состоЯНИЯ

= Spkzk + Bkпк. Данный алгоритм можно переформулировать в двухстадий-ной форме, т. е. разделить вычисления на два этапа: экстраполяции и фильтрации. Такое представление имеет особое значение в случае, когда измерения поступают через неравные промежутки времени. Такой вариант расширенного алгоритма обозначим как eSRCF-TM (eSRCF Time and Measurement Update).

В алгоритмах SRCF, SRCF-TM, SRIF, eSRCF и eSRCF-TM вычисления проводят с помощью матричного ортогонального преобразования QA = R оде A - в общем случае прямоугольная блочная матрица, Q - матрица QR-преобразования к верхней треугольной блочной матрице R Аналогично, в [20, р. 897] построен расширенный квадратно-корневой одностадийный информационный фильтр (алгоритм eSRIF).

eSRIF - extended Square-Root Information Filter

положить вычислить найти нач. оценку

Жо, Po = По > 0 Ro > 0 Qo > 0 det Fk = 0

Sn0, Sr0 , Sq0 - нижние треугольные факторы Холецкого

¿о = sno1xo

^^ Рекуррентное obhôbjiehpîe данных {к > 0^

-l Rk 0 0

-SR1Hk F,

l

S-lHkF— lGkSQk

-S—kZk

¿k 0

выполнить Qk Ak = Lk

Qk

s

k

1

>Pk

S -1F- 1

k l l -SpklFk lGkSQk

0

s

-1

Re,k — 1

Pk + i1 (*)

—S—1 Kp,k

s

0 1

k+i

(*)

0 0 W

i

-ëk

Zk+1 (*)

вычислить

sRk+i, sQ,

k+i

В нем в каждый момент времени k оценку вектора состояния вычисляют с помощью обратной подстановки в треугольной системе линейных алгебраических уравнений S—1 xk = zk. Избавиться от этих дополнительных вычислений можно путем добавления четвертого блочного столбца в алгоритм eSRIF [20, р. 898]. Такой алгоритм называют модифицированным квадратно-корневым информационным фильтром (алгоритм mSRIF - modified Square-Root Information Filter).

В алгоритмах eSRIF и mSRIF вычисления проводят с помощью матричного ортогонального преобразования QA = L, оде A - в общем случае прямоугольная блочная

матрица, 2 - матрица (^-преобразования к нижнему треугольному виду первых трех блочных столбцов матрицы А. Результатом преобразования является блочная матрица С.

2.1.2. Комбинированные ортогонализованные алгоритмы

Сочетание ковариационной и информационной форм в одном ортогонализован-ном алгоритме обосновано в работе [20]. Ниже приведены уравнения комбинированного квадратно-корневого фильтра сЭКР.

cSRF - combined Square-Root Filter

® Начальные данные

положить вычислить найти нач. оценку х0, Po = П0 > 0 (R > 0 Qo > 0 det Fk = 0) Sn0, Sr0, Sq0 - нижние треугольные факторы Холецкого Zo = Sp01xo

® Рекуррентное 0бн0вленр1е данных k > 0^ г oT n о-1гг тп-1 о-1гг тл-1^ а а-1.. п

Qk Ak = AÍ

Qk

Spk Hk Spk Fk 0 S~X, GT

SplFp 1 -SplFp 1<Gk SQk I

ST e,k 0 0 -e-k

0 S SPk+i Spk+i 01 Zk+1

0 0 (*) SQk (*)

Zk 0

вычислить

SRk+1, SQ

k+1

В алгоритме cSRF для блочной матрицы A вычисляют матричное ортогональное преобразование QA = A^ (QR-преобразование первых двух блочных столбцов либо QL-преобразование третьего и четвертого блочных столбца) и получают в результате блочную матрицу A

Комбинированный алгоритм cSRF является «универсальным> в том смысле, что позволяет одновременно вычислить все величины обоих алгоритмов - ковариационного SRCF и информационного SRIF. Вычислительная сложность алгоритма зависит от конкретной задачи и от тех величин фильтра, которые требуются для ее решения.

Алгоритм cSRF можно представить в двухстадийной форме в виде алгоритма cSRF-TM (Time and Measurement Update cSRF) [20, p. 899].

0

2.2. Алгоритмы на основе методов взвешенной ортогонализации

Второй класс ортогонализованных алгоритмов основан на методах взвешенной ортогонализации [21]. Известны два типа таких алгоритмов: UD-реализации и LD-реализации дискретного фильтра Калмана. Основными преимуществами алгоритмов являются устойчивость по отношению к ошибкам машинного округления, отсутствие операции извлечения квадратного корня, избавление от операции матричного обращения, компактная и однородная структура, ориентированность на параллельные вычисления.

Отличие данного класса алгоритмов заключается в том, что все они используют представление ковариационной матрицы Pk в виде матричного произведения CDCT, где C - верхняя треугольная U либо нижняя треугольная L матрица с единицами на D

ленной матрицы такое представление существует и может быть получено с помощью модифицированного разложения Холецкого [21]. Первой CI^модификацией дискретного фильтра Калмана был последовательный алгоритм Дж. Бирмана [13]8. Затем К. Торнтон [12] модифицировала метод взвешенной ортогонализации Грама - Шмидта и применила его для вычисления UD-факторов ковариационной матрицы Рк на этапе экстраполяции ФК.

Кратко изложим основную идею построения ортогонализованных алгоритмов на основе модифицированного метода взвешенной ортогонализации Грама - Шмидта (MWGS - Modified Weighted Gram-Shmidt orthogonalization). Рассмотрим матричное уравнение вида Р = ATDiA + CTD2C. Если существует матрица MWGS-B [ AT CT] = UB T , где U - верхняя треугольная мат-

рица с единицами на диагонали, то

= ATDiA + CtD2C = UDUT, (6)

где {U,D} верхний треугольный и диагональный факторы в UD-разложении матри-

Р основан на обратной процедуре взвешенной

ортогонализации при обновлении пары матриц {U,D} (6).

Рассмотрим современные ортогоналнзованные алгоритмы^ основанные на взве~ шенной ортогонализации.

2.2.1. Одностадийный ортогонализованный UD-алгоритм

В [16] авторы предложили ортогонализованный алгоритм этапа фильтрации дискретного ФК, основанный на быстрых вращениях Гивенса [21]. Предположим, что в результате выполнения этапа экстраполяции получены LD-факторы ковариационных матриц Рк = LPk Dpk L~T>k и Rk = LRk LRk L~Tk. Тогда на этапе фильтрации в алгоритме Калмана LD-факторы LPk]k и Dpk|k ковариационной матрицы Рк\к можно вычислить с помощью матричного ортогонального преобразования

LRk Нк LPk 0 LPk

где Qd - матрица быстрых вращений Гивенса, не содержащих операции извлечения квадратного корня, LRe k и DRek - LD-факторы ковариационной матрицы невязки измерений Re>к, Кк - матрица Калмана. В [16] показано, что данный алгоритм алгебраически эквивалентен LD-реализации алгоритма Бирмана для случая скалярных измерений [30, теорема 13.4].

В [23, с. 262] на основе указанных выше преобразований сформулирован одностадийный ортогонализованный UD-алгоритм. Вместо быстрых вращений Гивенса использован метод модифицированной взвешенной ортогонализации Грама - Шмидта. Назовем данный алгоритм UD-CF (UD Covariance Filter). Алгоритм UD-CF позволяет вычислить те же величины, что и алгоритм CKF, но в отличии от последнего, является устойчивым по отношению к ошибкам машинного округления. Вместе с тем, он обладает следующими недостатками: вычисление оценки вектора состояния

8Дж. Бирман не только доказал численную эффективность последовательного UD-фильтра, но и показал, что при соответствующей программной реализации его алгоритм не сложнее, чем обычный алгоритм Калмана, см. [11].

Р = [ AT CT ]

Di 0 A

0 D2 C

Q

D

L Re, k Кк LRe

0

LPk

(7)

kk

Хк происходит отдельно от обработки остальных данных, которые входят в блочные матрицы и обновляются процедурой М\\"(!Я-1>1> [огона, ш ;ашш. Таким образом, этот вариант алгоритма имеет неоднородную структуру.

Для преодоления указанных недостатков в недавней работе [34] построен расширенный ортогонализованный (Ю-алгоритм, в котором все данные обрабатываются единообразно с помощью процедуры М\\"(!Я-ор гогона. ш ;ашш.

2.2.2. Расширенный ковариационный ортогонализованный 1ГО-алгоритм

Идея расширенного ортогонализованного алгоритма заключается в том, что все величины дискретного фильтра помещаются в блочную матрицу и обрабатываются единообразно с помощью взвешенных ортогональных преобразований. Рассмотрим представление ковариационной матрицы в виде Рк = иРк БРк ир , где и -верхняя треугольная матрица, Б - диагональная матрица. Введем обозначения X = (иРкБРк)-1Хк, Ък = (иРе кк)-1вк. Также предположим, что в системе (1), (2) вектор управления ик = 0. Уравнения расширенного ортогонализованного 1ГО-фильтра запишем в следующем виде.

eUD-CF - extended UD Covariance Filter

® Начальные данные

нач. данные

вычислить

положить

Хо = Хо (По > 0 Ro > 0 Qo > 0) {Un0 ,Dn0 },{URo,DRo },{UQo ,DQo } Хо = (UnoDno) :Хо и UPo = Un0, D

UD-факторы Холецкого

Po = Dno

©Рекуррентное обновл ение ддтк (k > 0)

заполнить пред-массивы

ВЫПОЛНИТЬ MWGS

получить пост-массивы

ВЫЧИСЛИТЬ

AI

Dk AI

Ak "

Ak BI,

S

Rk+i >

t

1 0 0

SQk+i

0

Gk UQk 0

Diag{DQk

oT

Fk UPk

Hk UPk

-zI (URk DRk ) 0

URk

-T

DPk ,DRk}

где B1 Dk Bk

V

t

U

'k+l

Pk+1

bI

Kp,k Ur

URe

e, k

D\ = Diag{7, DPk+1 ,DRe, k}

Алгоритм е1ГО-СР легко реализовать в виде программного кода на ЭВМ, поскольку он имеет компактную и однородную структуру, позволяющую единообразно обрабатывать все данные. В алгоритме отсутствует операция извлечения квадратного корня, а вычисление оценки вектора состояния хк не требует операции матричного обращения, которая, как известно, является чувствительной к ошибкам машинного округления.

0

2.2.3. Другие варианты алгоритмов

Объединяя алгоритм Торнтон для ниж^них треугольных факторов [30, с. 279] с блочным представлением Кайлата (7), в [38] впервые представлен двухстадииныи ор~ тогонализованного алгоритм LD-TM (LD Time and Measurement Update). Здесь на этапе экстраполяции применяют прямую процедуру MWGS-ортогонализации строк матрицы W с весовой матрицей D и получают в результате нижнюю треугольную

блочную матрицу ^^ и блочно-диагональную матрицу ; вычисления на этапе фильтрации проводят аналогично.

В данном разделе рассмотрены только алгоритмы ковариационного типа. Что касается алгоритмов информационного типа, то авторам известна только одна работа [39], в которой построен один из возможных вариантов информационного 1ГО-фильтра. Однако, этот алгоритм не обладает удобным для реализации блочно-матричным представлением как, например, алгоритм е1ГО-СЕ, поскольку в нем все необходимые 1ГО-факторы ковариационных матриц вычисляются по отдельности с помощью М\¥С8-преобразования и с помощью метода одноранговой модификации. Оценки вектора состояния также вычисляются по отдельным выражениям.

Таким образом, в настоящее время актуальной является задача построения новых ортогонализованных 1ГО-фильтров информационного типа и исследование их вычислительных свойств в сравнении с существующими информационными реализациями ФК.

2.3. 3-ортогонализованные быстрые алгоритмы

Третий важный класс ортогонализованных алгоритмов построен на базе 3-ортогональных

преобразований. Е го называют классом быстрых алгоритмов Чанд-расекара9 [8-10,19]. Их достоинствами являются устойчивость по отношению к ошибкам округления, свойственная всем ортогонализованным блочным алгоритмам, ориентированность на параллельные вычисления, существенное сокращение вычислительной сложности.

Поскольку применение ортогональных преобразований позволяет строить численно устойчивые модификации фильтра Калмана, ортогонализованные алгоритмы были получены и для класса

быстрых алгоритмов. Однако, необходи-

мо обновлять уравнения вида СТС = Лт Л ± БтБу что невозможно осуществить с помощью обычных ортогональных преобразований. С этой целью используют 3-ортогональные преобразования вида QЛ = Я, где матрица преобразования Q является 3-ортогональной [25].

Рассмотрим ортогонализованную квадратно-корневую форму быстрых алгоритмов [22].

J-SRFF - Square-Root J-orthogonal Fast Filter

нач. значения вычислить найти разложение Жс|о = жо, Рою = По, По > 0 Кро = ^ПоЯтЕ"о1/2, Re,о = R + ЯПоЯт 0Ро = ЬоЗЬт, где Lо € rraxa, S € raxa - сигнатура матрицы 0Ро и а = rank(5Pb), 5Ро = F^Fт + GQGT - КР,оК!о - По

® Обновление матриц фильтр а (к > 0)

вычислить R1/2 KT R1/2 KT Q, Re,k Kp,k = Re,k+1 Kp,k+1 Qk LTHT L TF t 0 L T+1

® Фильтрация (к > 0)

оценка ковариация ошибки Xk+1 = Fxk + Buk + Kp,k ek k Pk+1 = По + £ L j LT з=о

9От англ.: fast Chandrasekhar-type algorithms. Полное название таких методов Chandrasekhar-Kailath-Morf-Sidhu алгоритмы.

В общем случае, матрицы R не вырождены И Н6 обязательно ЯВЛЯЮТСЯ ЗНЭ/~ неопределенными. Однако они имеют очень важное для построения ортогонализо-ванных квадратно-корневых алгоритмов свойство, а именно, постоянную во времени энерцию. Рассмотрим разложение вида Rr,k = -RJ/2SkRjk-> àPk = Pk+i\k — Pk\k-l =

LkSkLJ, где Lk = LkRr]J2, Lk E Rraxa и Sk E Raxa - сигнатура эрмитовой матрицы, т.е. диагональная матрица, в которой количество ±1 на главной диагонали определено числом положительных и отрицательных собственных значений матрицы SPk-В силу закона энерции эрмитовых форм Sk = S, Ук и S - сигнатура матрицы ôPq. В алгоритме SRFF Qk - матрица J-ортогонального преобразования (сигнатура J = Im ф S = diag{Im,S}), приводящего к верхнему треугольному виду первый блочный столбец матрицы, стоящей в левой части матричного уравнения шага ® в алгоритме SRFF, т. е. RJk2+l - верхняя треугольная матрица. Теоретическое обоснование данного класса алгоритмов см. в [22, гл. 13].

2.4. Алгоритмы на основе сингулярного разложения

Четвертым классом современных ортогонализованных алгоритмов являются алгоритмы, построенные на основе сингулярного разложения (SVD-преобразования). Главная особенность заключается в том, что для их реализации не требуется строгая

Qk

Rk, а также начальной матрицы ковариации о шибки оценки P0. Другими словами, если в процессе работы такого алгоритма положительная определенность указанных матриц будет потеряна, это обстоятельство никак не повлияет на работоспособность SVD-фильтра.

Идея построения SVD-фильтра состоит в представлении ковариационных матриц ошибок оценивания в виде Pk = Opfc Dpk Opfc, где Opk - ортогональная матрица и Dpk -

Pk

Pk

щью процедуры SVD-факторизации (сингулярного разложения). Насколько известно авторам данной статьи, впервые SVD-фильтр был построен в [18, р. 1226]. Кроме сингулярного разложения, в уравнениях фильтра использовалось также разложение Холецкого для вычисления квадратно-корневых факторов матриц ковариаций шу-Qk Rk

вычисляют сингулярное разложение прямоугольной блочной матрицы Л, в результате получают SVD-факторы ОТ, S и V такие, что Л = OTSVт. Следует отметить, что алгоритм SVD-SRKF не свободен от недостатков. Во -первых, на этапе фильтрации требуется три операции матричного обращения. Во-вторых, реализация алгоритма

Qk Rk

причинам (условия задачи или вычисленные в процессе работы алгоритма значения) эти матрицы потеряют положительную определенность, разложение Холецкого невозможно будет выполнить, и алгоритм SVD-SRKF потеряет работоспособность.

Результатом устранения указанных недостатков стал улучшенный SVD-фильтр, предложенный в недавней работе авторов [36]. Уравнения алгоритма представлены ниже.

SVD-KF - improved SVD Kaiman Filter

вычислить НУ!) По = 0по Dno 0Ro q0 = 0Qo DQo 0Qo

положить x0|0 = x0, 0Po\o = 0по, Dl/2 = D^2 Po\o no

© Экстраполяция (к > 0)

D/2e l.. FR

вычислить НУ!)

найти оценку вычислить НУ!)

P.

k\k

DQk 0Qk Gк

= Wtu

xk+1 = Fk xk\k + Bk uk Qk+1 = 0Qk+1 DQk+1 &Qk+1,

D

1/2 Pk+1 0

0jk+1

Rk+1 = 0R,

k+1

DRk+10jk+1

©Фильтрация (k > 0)

DRk + 1 0 jk + 1 DPk+1 0 Rk+1 Hk+H

вычислить НУ!)

наити

вычислить НУ!)

наити оценку

W

(1) MU

D

1 2

Re

,k+1

0R

Re,k + 1

Kk+1 = ifc+1hj+i0Re,k+1,

D^L 0 1+1 (i " kk+hk+1)

Kk+1=Kk+1DR:,k+10 ik+1

T■

DRi+10 ~Rk+1Kk+1

W

(2) MU

-1 1/2

D

Pk

k+1\k+1 0

0

T

Pk

k+1|k+1

xk+1\k+1 = xk+1 + Kk+1d—1 k+1 ek+1

ek+1 = 0Rek+1 (zk+1 - Hk+1xk+1)

В отличие от ЯУ1)-ЯИКК. в алгоритме ЯУ1)-КК требуется всего одно обращение диагональной матрицы О—, . В [36] доказано, что алгоритм ЯУБ-КР алгебраически эквивалентен стандартному представлению фильтра Калмана СКК-ТМ. НО В ОТЛИчие от последнего, обладает улучшенными вычислительными свойствами, а именно, высокой устойчивостью к ошибкам машинного округления. Алгоритм ЯУ1)-КК отличается от других ортогонализованных алгоритмов тем, что в отличие от ЯК- и 1ГО-алгоритмов он не требует условия положительной определенности матриц кова-риаций шумов Qk и и в отличие от алгоритмов информационного типсЦ не требует невырожденности переходной матрицы состояния Fk. Следовательно, данный алгоритм может применяться в качестве численно эффективной реализации фильтра Калмана для широкого класса задач дискретной фильтрации.

Также, как и в предыдущем пункте, здесь рассмотрены ЯУБ-реализацпп ковариационного типа. Насколько известно авторам данной статьи, для алгоритмов информационного типа ЯУБ-реалпзацпй пока не существует, и эта задача может быть предметом отдельного исследования.

0

3. Вычислительные аспекты ортогонализованных алгоритмов

Рассмотрим вычислительные аспекты современных ортогонализованных алгоритмов. Традиционно, устойчивые реализации дискретного фильтра Калмана подразделяют на следующие категории: алгоритмы ковариационного и информационного типа (обновляющие на каждом шаге, соответственно, значения ковариационной матрицы Р либо значения информационной матрицы Л); двухстадииные и одностадийные (совмещенные) алгоритмы; расширенные алгоритмы.

В [37] предложена классификация современных ортогонализованных алгоритмов в зависимости от типа матричных ортогональных преобразований (см. таблицу).

Основным критерием классификации является способ построения и численной реализации ортогонализованного алгоритма. Приведенный обзор показал, что каж-

Таблица

Классификация современных ортогонализованных алгоритмов

Ортогонализованный алгоритм Тип ортогонального преобразования

Алгоритмы ковариационного типа

® SRCF [5-7] Q Ai A2 = R 0

© SRCF-TM [5-7] qTU A A 11 12 = R1 [0 ; ^UA2 = R2

® eSRCF [20,22] Q A A ii Ai2 2i A22 = R11 R12 [0 R22

® eSRCF-TM [20,22] QTUAi = Ri; QMU A2i A22 — R21 R22

® UD-CF [12,13,16] A 1 = UB 1 , где A 1 Dw A — UDß U 1

© eUD-CF [23,34] A 1 — UB 1 , где A 1 Dw A — UDß U 1

© LD-TM [30,38] ATU = LiBT , где (ATU) 1 D^ ATU = LxDßlL1 ; AMU = L2Bj, где (AMU)TDWVAmu = L2Dß2LT2

® J-SRFF [9,10,19,22] QA — R, где Q 1 JQ — J

® SVD-SRKF [18] A'rи = wti . amu = 1

® SVD-KF [36] дги = ^T 1 . Afv = SV t', Afv = W2MU SV T

Алгоритмы информационного типа

® SRIF [4,13] Qi4J Aii Ai2 = R11 R12 ; ^UA2 — R2

© eSRIF [20,22] Q A21 A22 = [ L2i L22 i

® mSRIF [20,22] Q A2i A22 = L2i L22

Комбинированные алгоритмы

® cSRF [20] QA = [R|L]

© cSRF-TM [20] ^UA = [R|L^ ^MUA = [R|L]

Всего: 15 классов алгоритмов

дый из алгоритмов основан на применении некоторого матричного ортогонального преобразования и не зависит от способа его вычисления.

Таким образом, выбирая различные методы ортогональных преобразований, можно построить целые семейства ортогонализованных алгоритмов. Более того, зная принцип построения ортогонализованного алгоритма на основе того или иного матричного ортогонального преобразования, можно построить и теоретически обосновать новые ортогонализованные алгоритмы ковариационного и информационного типов, а также ортогонализованные алгоритмы для других классов дискретных систем, например, для дискретных линейных систем с коррелированными шумами, для нелинейных дискретных систем, для парных марковских моделей систем и других.

3.1. Особенности программной реализации ортогонализованных алгоритмов

В таблице показано, что каждый ортогонализованный алгоритм построен на основе одного из матричных ортогональных преобразований. Следовательно, методы ортогональных преобразований можно рассматривать как базовые вычислительные технологии для реализации ортогонализованных алгоритмов. Поэтому технология программной реализации ортогонализованных алгоритмов представляется простой и

понятной [37]. Она состоит из трех основных этапов:

1. Заполнить имеющимися данными блочньтй(ьте) пред-массив(ьт) [ргес1-аггау(з)].

2. Выполнить матричное ортогональное преобразование одного из типов, представленных в табл.

3. Получить в результате блочньтй(ые) пост-массив (ьт) [ров^аггау^)] и извлечь из него (из них) в виде матричных блоков искомые результаты.

Общие схемы программной реализации различных классов алгоритмов показаны на рис. 1 и 2 а). Примеры использования этих общих схем показаны на рис. 1 и 2 б).

а) б)

Рис. 1. а) Общая схема программной реализации квадратно-корневых алгоритмов; б) Пример схемы программной реализации алгоритма БШР (этап фильтрации) па основе ОБ, преобразования

а) б)

Рис. 2. а) Общая схема программной реализации СБ алгоритмов; б) Пример схемы программной реализации алгоритма 1Л)-СР па основе МД¥СБ-иВ преобразования

Аналогично можно представить общие схе м ы дл я С^Ь-. ЬБ-. ,ТС^К- и БУБ-преобразований. Заполняя блочные массивы необходимыми данными, получаем конкретный вариант ортогонализованного алгоритма. Рис. 3 содержит пример программного кода этапа фильтрации алгоритма БЫЕ. Рис. 4 содержит пример программного кода алгоритма СБ-СЕ.

Для реализации метода (^Я-преобразования использована т-функция qr из стандартной библиотеки МАТЬАВ и программная реализация в виде т-функции

function [R_L,zx] = SRIF_TM_update(Н,S_R,z,R_L,zx) [m,n] = size(H); A = [R_L zx; S_R\H S_R\z]; [Q,RR] = qr(A); R_L = RR(1:n,l:n); zx = RR(l:n,n+l); end

Рис. 3. Программный код этапа фильтрации алгоритма SRIF

function [x,U_P,D_P] = UD_CF(F,В>G>Н,U_Q,D_Q,U_R,D_R,U_P,D_P,z,x,u) [m,n] = size(H); q = size(G,2);

A = [G*U_Q F*U_P zeros(n,m); zeros(m,q) H*U_P U_R]';

Dw = diag([diag(D_Q); diag(D_P); diag(D_R)]);

[U, Db] = MWGS_UD(A, Dw);

U_P = U(l:n,l:n);

D_P = Db(l:n,l:n);

\

end

Рис. 4. Программный код алгоритма Г1 )-(Т"

МИСЗ.ГО обратной процедуры модифицированной взвешенной ортогонализации Грама - Шмидта [13].

Таким образом, программная реализация ортогонализованных алгоритмов представляет собой применение матричных ортогональных преобразований к блочным массивам данных. Такая однородная структура алгоритмов идеально подходит для распараллеливания вычислений и реализации средствами параллельного программирования.

3.2. Численная устойчивость алгоритмов к ошибкам машинного округления

Проблема влияния ошибок машинного округления на работоспособность стандартного алгоритма Калмана была обозначена выше в разделе 1.2. Покажем на одном практическом примере улучшенные вычислительные свойства современных ортогонализованных алгоритмов по сравнению со стандартным алгоритмом Калмана.

Пример 1. Рассмотрим задачу оценивания параметров спутника, движущегося по круговой орбите [3, р. 1448]. Модель движения описывается разностными уравнениями

Хк+1

1 1 0, 5 0, 5 0 0 0 0

0 1 1 1 xk+wk, Q = 0 0 0 0

0 0 1 0 0 0 0 0

0 0 0 0, 606 0 0 0 q

, q = 0, 63 • 10

-2

(8)

с начальными условиями x0 ~ N(0, П0), П0 = I4.

Чтобы показать отличие в вычислительных свойствах алгоритмов дискретной фильтрации, рассмотрим плохо обусловленную схему измерений [23, р. 205]:

гк

111 1 1 1 1 1+8

Хк + Ук, Ук (0,Я),

где Я = 82/2. Для моделирования ситуации появления ошибок машинного округления предположим, что 82 < е0кр> но 8 > е0кр> гДе £окр _ параметр машинного округления.

Для иллюстрации потери точности вычислений стандартным фильтром Калма-на в условиях примера 1, приведем здесь результаты сравнительного исследования, проведенного авторами данной статьи и представленного в недавней работе [36].

Вычислительные эксперименты выполнены в системе матьав, обеспечивающим представление вещественных чисел с погрешностью 10_16. Параметр машинного округления еокр хранится в переменной ерэ, со значением еокр=ерз/2 = 1,1102 • 10_16. Мы провели серию вычислительных экспериментов для различных значений параметра 8 таких, что 8 ^ £ощ>. Методика проведения экспериментов следующая. Начиная с начального значения, для к = 1,... ,К = 100 мы провели компьютерное моделирование «истинной> траектории движения объекта в соответствии с (8). В результате были получены наборы значений вектора состояния ХкХаС и соответствующих измерений гк, к = 1,...,К. Далее, для сравнения различных алгоритмов, с помощью т-функций на языке Матьав мы посчитали значения оценок Хщ*, к = 1,..., К вектора состояния системы. Для всех рассмотренных алгоритмов фильтрации условия проведения эксперимента были одинаковые.

В целом, мы провели серию из М = 500 экспериментов. При каждом запуске вычислялась «истинная> траектория движения объекта, хк*ас1:, к = 1,... ,К, моделировались соответствующие измерения гк, к =1,...,Ки вычислялись значения оценок вектора состояния Хк\к, к = 1,...,К. По результатам экспериментов для каждого компонента вектора состояния мы вычислили среднеквадратичную ошибку (И МЯК) по формуле:

ЕМЯЕа

\

1 м к 2

^ (Х'к - Х^к\к) , (9) 3 = 1 к=1

где М = 500 - число запусков экспериментов, К = 100 - число измерений, х3'кхас и Х3к\к г-й компонент «истинного> (модельного) вектора состояния и его оценки, полученной при проведении j-тo эксперимента. Результаты экспериментов представлены на рис. 5.

По рис. 5 видно, что стандартный алгоритм Калмана (СКК) быстро теряет точность вычислений при 8 ^ £0кр- При 8 = 10_7 и далее, стандартный метод уже не пригоден для вычислений, поскольку ||КМ8ЕХ;||2 = КаК, т.е. из-за ошибок машинного округления невозможно получить результат (^для наглядности мы отобразили значения МаМ на верхней границе графика). Напротив, ортогонализованным алгоритмам удается поддерживать приемлемую точность вычислений вплоть до значения 8 = 10"14. Подробный сравнительный анализ рассмотренных алгоритмов см. в [36].

Таким образом, в примере 1 показано, что современные ортогонализованные алгоритмы являются наименее чувствительными к ошибкам машинного округления, и следовательно, наиболее надежными алгоритмами оптимального линейного оценива-

Заключение

В настоящее время вычислительные методы оптимального оценивания стали самостоятельной областью исследования и получили большое развитие. Современные численно эффективные ортогонализованньте алгоритмы привлекательны не только своей устойчивостью к ошибкам машинного округления. Важным качеством, которое отмечается в научной литературе, является приспособленность алгоритмов, использующих различные типы матричных ортогональных преобразований, к программной реализации на параллельных или векторных вычислительных системах, что позволяет разрабатывать новые эффективные информационные технологии, в частности, при решении задач реального времени и при обработке больших дан Н ЫХ •

Любая ортогонализованная реализация дискретного фильтра Калмана (как в стандартной5 так и в быстрой форме) предполагает использование ортогональных преобразований на каждом тттаге алгоритма. Метод ортогонализации и. следовательно. матрицы ортогонального преобразования допускают определенный произвол, т. е. для реализации того или иного алгоритма можно использовать любые ортогональные преобразования, приводящие к желаемому виду блочной матрицы. Следовательно, вычислительная сложность и численная эффективность ортогонализованньтх алгоритмов определяется в основном сложностью и вычислительными свойствами применяемых методов ортогонализации.

В данной статье дан краткий обзор существующих классов современных ортогонализованньтх алгоритмов оптимальной дискретной фильтрации. Предложенная классификация алгоритмов по типам матричных ортогональных преобразований позволяет легче понять метод их построения и выбрать способ эффективной программной реализации при решении практических задач оценивания в классе многомерных дискретных линейных стохастических систем.

Все рассмотренные алгоритмы являются алгебраически эквивалентными стандартным реализациям дискретного фильтра Калмана, но существенно превосходят его по своим вычислительным свойствам. Поэтому применение ортогонализованных алгоритмов при решении практических задач позволяет получить численно эффективные и надежные решения.

М.В. Куликова благодарит португальский Фонд науки и технологии (Fundagao para a Ciencia е a Tecnología) за оказанную финансовую поддержку в рамках проекта

UID/МиШ/04 621/2013. Литература

1. Kalman, R.E. A New Approach to Linear Filtering and Prediction Problems / R.E. Kalman // Journal of Basic Engineering. - 1960. - V. 82, № 1. - P. 35-45.

2. Potter, J.E. Statistical Filtering of Space Navigation Measurements / J.E. Potter, R.G. Stern // Proceedings 1963 AIAA Guidance and Control Conference. - 1963. - 13 p.

3. Rauch, H.E. Maximum Likelihood Estimates of Linear Dynamic Systems / H.E. Rauch, C.T. Striebel, F. Tung // AIAA journal. - 1965. - V. 3, № 8. - P. 1445-1450.

4. Dyer, P. Extension of Square-Root Filtering to Include Process Noise / P. Dyer, S. McReynolds // Journal of Optimization Theory and Applications. - 1969. - № 3. -P. 444-459.

5. Schmidt, S.F. Computational Techniques in Kalman Filtering / S.F. Schmidt // Theory and Applications of Kalman Filtering. NATO Advisory Group for Aerospace Research and Development, AGARDograph 139. - 1970.

6. Kaminski, P.G. Discrete Square-Root Filtering: A Survey of Current Techniques / P.G. Kaminski, A.E. Bryson, S.F. Schmidt // IEEE Transactions on Automatic Control. -1971. - V. 16, № 6. - P. 727-735.

7. Kaminski, P.G. Square Root Filtering and Smoothing for Discrete Processes. Ph.D. Thesis / P.G. Kaminski. - Stanford University, 1971.

8. Kailath, T. Some New Algorithms for Recursive Estimation in Constant Linear System / T. Kailath // IEEE Transactions on Information Theory. - 1973. - V. 19, № 11. - P. 750-760.

9. Morf, M. Some New Algorithms for Recursive Estimation in Constant, Linear Discrete-Time Systems / M. Morf, G.S. Sidhu, T. Kailath // IEEE Transactions on Automatic Control. -1974. - V. 19, № 4. - P. 315-323.

10. Morf, M. Square-Root Algorithms for Least-Squares Estimation /М. Morf, T. Kailath // IEEE Transactions on Automatic Control. - 1975. - V. 20, № 4. -P. 487-497.

11. Thornton, C.L. Numerical Comparison of Discrete Kalman Filtering Algorithms: an Orbit Determination Case Study / C.L. Thornton, G.J. Bierman // JPL Technical Memorandum 33-771. - Pasadena, 1976. - 48 p.

12. Thornton, C.L. Triangular Covariance Factorizations for Kalman Filtering. Ph.D. Thesis / C.L. Thornton. - University of California at Los Angeles, 1976.

13. Bierman, G.J. Factorization Methods For Discrete Sequential Estimation / G.J. Bierman. -New York: Academic Press, 1977.

14. Фомин, B.H. Рекуррентное оценивание и адаптивная фильтрация / В.Н. Фомин. - М.: Наука, 1984.

15. Verhaegen, М. Numerical Aspects of Different Kalman Filter Implementations / M. Verhaegen, P. Van Dooren // IEEE Transactions on Automatic Control. - 1986. - V. 31, № 10. - P. 907-917.

16. Jover, J.M. A Parallel Architecture for Kalman Filter Measurement Update and Parameter Estimation / J.M. Jover, T. Kailath // Automatica. - 1986. - V. 22, № 1. - P. 43-57.

17. Огарков, M.A. Методы статистического оценивания параметров случайных процессов / М.А. Огарков. - М.: Энергоатомиздат, 1990.

18. Wang, L. Kalman Filter Algorithm Based on Singular Value Decomposition / L. Wang, G. Libert, P. Manneback // Proceedings of the IEEE Conference on Decision and Control, Tuczon, 1992. - P. 1224-1229.

19. Sayed, A.H. Extended Chandrasekhar Recursions / A.H. Sayed, T. Kailath // IEEE Transactions on Automatic Control. - 1994. - V. 9, № 3. - P. 619-622.

20. Park, P. New Square-Root Algorithms for Kalman Filtering / P. Park, T. Kailath // IEEE Transactions on Automatic Control - 1995. - V. 40, № 5. - P. 895-899.

21. Голуб, Дж. Матричные вычисления / Дж. Голуб, Ч. Ван Доун. - М.: Мир, 1999.

22. Kailath, Т. Linear Estimation / Т. Kailath, A.H. Sayed, В. Hassibi. - New Jersey: Prentice Hall, 2000.

23. Grewal, M.S. Kalman Filtering: Theory and Practice Using MATLAB / M.S. Grewal, A.P. Andrews. - New Jersey: Prentice Hall, 2001.

24. Higham, N.J. Accuracy and Stability of Numerical Algorithms / N.J. Higham. - Philadelphia: SI AM, 2002.

25. Higham, N.J. J-Orthogonal Matrices: Properties and Generalization / N.J. Higham // SIAM Review. - 2003. - V. 45, № 3. - P. 504-519.

26. Zhou, B.A. UB Factorization-Based Nonlinear Adaptive Set-Membership Filter for Ellipsoidal Estimation / B. Zhou, J. Han, G. Liu // International Journal of Robust and Nonlinear Control. - 2008. - V. 18. - P. 1513-1531.

27. Baowang, F. Square-Root Second-Order Extended Kalman Filter and its Application in Target Motion Analysis / F. Baowang, L. Teng, H.Z. Tao // Radar, Sonar and Navigation. - 2010. - V. 4, № 3. - P. 329-335.

28. Gibbs, B.P. Advanced Kalman Filtering, Least-Squares and Modeling: A Practical Handbook / B.P. Gibbs. - Hoboken, New Jersey: John WTiley and Sons, 2011.

29. Семушпн, И.В. Адаптивные системы фильтрации, управления и обнаружения: коллективная монография / И.В. Семушпн, Ю.В. Цыганова, М.В. Куликова и др. - Ульяновск: УлГУ, 2011.

30. Семушпн, И.В. Вычислительные методы алгебры и оценивания / И.В. Семушин. - Ульяновск: УлГТУ, 2011.

31. Семушин, И.В. Устойчивые алгоритмы фильтрации - обзор и новые результаты для систем судовождения / И.В. Семушин, Ю.В. Цыганова, К.В. Захаров // Информационные технологии и вычислительные системы. - 2013. - № 4. - С. 90-112.

32. Kulikova, M.V. Stochastic Volatility Models for Exchange Rates and Their Estimation Using Quasi-Maximum-Likelihood Methods: an Application to the South African Rand / M.V. Kulikova, B.R. Taylor // Journal of Applied Statistics. - 2013. - V. 40, № 3. -P. 495-507.

33. Kulikova, M.V. Adaptive Wave Filtering for Marine Vessels within UB-Based Algorithms / M.V. Kulikova, J.V. Tsyganova, I.V. Semushin // Proceedings of the ECC2016, European Control Conference. - Aalborg, 2016. - P. 807-812.

34. Цыганова, Ю.В. О методах реализации UD-фильтра / Ю.В. Цыганова // Известия высших учебных заведений. Поволжский регион. Физико-математические науки. - 2013. -№ 3. - Р. 84-104.

35. Semushin, I.V. Numerically Efficient UD Filter Based Channel Estimation for OFDM Wireless Communication Technology / I.V. Semushin, Yu.V. Tsyganova, A.V. Tsyganov, E.F. Prokhorova // Procedia Engineering. - 2017. - V. 201. - P. 726-735.

36. Kulikova, M.V. Improved Discrete-Time Kalman Filtering within Singular Value Decomposition / M.V. Kulikova, J.V. Tsyganova // IET Control Theory and Applications. - 2017. - V. 11, № 5. - P. 2412-2418.

37. Цыганова, Ю.В. Ортогонализованные блочные методы для параметрической идентификации дискретных линейных стохастических систем: дис. ... д-ра физ.-мат. наук / Ю.В. Цыганова. - Ульяновск, 2017.

38. Semushin, I.V. Numerically Efficient LB-computations for the Auxiliary Performance Index Based Control Optimization under Uncertainties / I.V. Semushin, J.V. Tsyganova, A.V. Tsyganov // Proceedings of the 17th IFAC Workshop on Control Applications of Optimization, Yekaterinburg, Russia, October 15-19, 2018. - Yekaterinburg, 2018. - V. 51. -P. 568-573.

39. B'Souza, Ch. Information Formulation of the UBU Kalman Filter / Ch. B'Souza, R. Zanetti // IEEE Transactions on Aerospace and Electronic Systems. - 2018.

Юлия Владимировна Цыганова, доктор физико-математических наук, профессор, кафедра «Информационные технологии:», Ульяновский государственный университет (г. Ульяновск, Российская Федерация), [email protected].

Мария Вячеславовна Куликова, К^НДИДсХТ физико-математических наук, научный сотрудник, СЕМАТ, Высший технический институт, Университет Лиссабона (г. Лиссабон, Португалия), [email protected].

Поступила в редакцию 17 августа 2018 г.

MSC 93A30, 93E10 DOI: 10.14529/mmp180401

ON MODERN ARRAY ALGORITHMS FOR OPTIMAL DISCRETE FILTERING

Yu.V. Tsyganova1, M.V. Kulikova2

Ulyanovsk State University, Ulyanovsk, Russian Federation 2Instituto Superior Tecnico, Universidade de Lisboa, Lisboa, Portugal E-mails: [email protected], [email protected]

Nowadays, computational methods for optimal estimation have become an independent field of research and have received a great progress. Modern numerically efficient array algorithms are attractive not only because of their robustness to machine round-off errors, but additionally because of utilization various types of matrix orthogonal transformations. Thus, their design pattern is well suited for parallel implementations on modern computing systems. These properties allow to develop new efficient information technologies, in particular, the techniques that are applicable for solving real-time problems as well as for processing big data arrays. This paper gives a brief survey of modern array algorithms for optimal linear discrete-time filtering. Four large classes of array algorithms are considered:

square-root array algorithms, array algorithms based on weighted orthogonalization, J-orthogonal array algorithms and methods based on singular value decomposition. We suggest a classification of array algorithms according to the type of the utilized matrix orthogonal transformation on the basis of which these algorithms are designed. Such classification suggests a more simple way for understanding the array filtering methods' design and gives a choice for finding their most effective implementation for estimating multivariable discrete-time linear stochastic systems. The computational aspects of array algorithms are investigated. It includes the numerical stability to machine round-off errors, and discussion of efficient software implementation for the array algorithms under examination. Finally, the array algorithms investigated in this paper are algebraically equivalent to the conventional implementation of the discrete-time Kalman filter, but they possess the significantly improved computational properties. The results of the presented comparative study allow to conclude that the use of array algorithms in solving practical problems helps to obtain numerically efficient and reliable solutions.

Keywords: discrete filtering; linear stochastic systems; Kalman filter; matrix orthogonal transforms; array algorithms.

References

1. Kalman R.E. A New Approach to Linear Filtering and Prediction Problems. Journal of Basic Engineering, 1960, vol. 82, no. 1, pp. 35-45. DOI: 10.1115/1.3662552

2. Potter J.E., Stern R.G. Statistical Filtering of Space Navigation Measurements. AIAA Guidance, Navigation and Control Conference, Cambridge, 1963, 13 p. DOI: 10.2514/6.1963333

3. Rauch H.E., Striebel C.T., Tung F. Maximum Likelihood Estimates of Linear Dynamic Systems. AIAA journal, 1965, vol. 3, no. 8, pp. 1445-1450. DOI: 10.2514/3.3166

4. Dyer P., McReynolds S. Extension of Square-Root Filtering to Include Process Noise. Journal of Optimization Theory and Applications, 1969, no. 3, pp. 444-459. DOI: 10.1007/BF00929358

5. Schmidt S.F. Computational Techniques in Kalman Filtering Theory and Applications of Kalman Filtering. NATO Advisory Group for Aerospace Research and Development, 1970.

6. Kaminski P.G., Bryson A.E., Schmidt S.F. Discrete Square-Root Filtering: A Survey of Current Techniques. IEEE Transactions on Automatic Control, 1971, vol. 16, no. 6, pp. 727-735. DOI: 10.1109/TAC.1971.1099816

7. Kaminski P.G. Square Root Filtering and Smoothing for Discrete Processes. Ph.D. Thesis, Stanford University, 1971.

8. Kailath T. Some New Algorithms for Recursive Estimation in Constant Linear System. IEEE Transactions on Information Theory, 1973, vol. 19, no. 11, pp. 750-760. DOI: 10.1109/TIT. 1973.1055104

9. Morf M., Sidhu G.S., Kailath T. Some New Algorithms for Recursive Estimation in Constant. Linear Discrete-Time Systems. IEEE Transactions on Automatic Control, 1974, vol. 19, no. 4, pp. 315-323. DOI: 10.1109/TAC.1974.1100576

10. Morf M., Kailath T. Square-Root Algorithms for Least-Squares Estimation. IEEE Transactions on Automatic Control, 1975, vol. 20, no. 4, pp. 487-497. DOI: 10.1109/TAC. 1975.1100994

11. Thornton C.L., Bierman G.J. Numerical Comparison of Discrete Kalman Filtering Algorithms: an Orbit Determination Case Study. JPL Technical Memorandum 33-771, 1976, 48 p.

12. Thornton C.L. Triangular Covariance Factorizations for Kalman Filtering. Ph.D. Thesis, University of California at Los Angeles, 1976.

13. Bierman G.J. Factorization Methods For Discrete Sequential Estimation. N.Y., Academic Press, 1977.

14. Fomin V.N. Rekurrentnoe otsenivanie i adaptivnaya fil'tratsiya [Recurrent Estimation and Adaptive Filtering]. Moscow, Nauka, 1984. (in Russian)

15. Verhaegen M., Van Booren P. Numerical Aspects of Bifferent Kalman Filter Implementations. IEEE Transactions on Automatic Control, 1986, vol. 31, no. 10, pp. 907-917. BOI: 10.1109/TAC.1986.1104128

16. Jover J.M., Kailath T. A Parallel Architecture for Kalman Filter Measurement Update and Parameter Estimation. Automatica, 1986, vol. 22, no. 1, pp. 43-57. BOI: 10.1016/0005-1098(86)90104-4

17. Ogarkov M.A. Metody statisticheskogo otsenivaniya parametrov sluchajnykh protsessov [Methods of Statistical Estimation of Parameters of Random Processes]. Moscow, Energoatomizdat, 1990. (in Russian)

18. Wang L., Libert G., Manneback P. Kalman Filter Algorithm Based on Singular Value Becomposition. Proceedings on the IEEE Conference on Decision and Control, 1992, pp. 1224-1229. BOI: 10.1109/CBC.1992.371522

19. Sayed A.H., Kailath T. Extended Chandrasekhar Recursions. IEEE Transactions on Automatic Control, 1994, vol. 9, no. 3, pp. 619-622. BOI: 10.1109/9.280773

20. Park P., Kailath T. New Square-Root Algorithms for Kalman Filtering. IEEE Transactions on Automatic Control, 1995, vol. 40, no. 5, pp. 895-899. BOI: 10.1109/9.384225

21. Golub G., Van Loan Ch. Matrix Computations. Johns Hopkins University Press, 1996.

22. Kailath T., Sayed A.H., Hassibi B. Linear Estimation. New Jersey, Prentice Hall, 2000.

23. Grewal M.S., Andrews A.P. Kalman Filtering: Theory and Practice Using MATLAB. New Jersey, Prentice Hall, 2001.

24. Higham N.J. Accuracy and Stability of Numerical Algorithms. Philadelphia, SIAM, 2002.

25. Higham N.J. J-Orthogonal Matrices: Properties and Generalization. SIAM Review, 2003, vol. 45, no. 3, pp. 504-519. BOI: 10.1137/S0036144502414930

26. Zhou B., Han J., Liu G. A UB Factorization-Based Nonlinear Adaptive Set-Membership Filter for Ellipsoidal Estimation. International Journal of Robust and Nonlinear Control, 2008, vol. 18, pp. 1513-1531. BOI: 10.1002/rnc.l289

27. Baowang F., Teng L., Tao H.Z. Square-Root Second-Order Extended Kalman Filter and Its Application in Target Motion Analysis. Radar, Sonar and Navigation, 2010, vol. 4, no. 3, pp. 329-335. BOI: 10.1049/iet-rsn.2008.0070

28. Gibbs B.P. Advanced Kalman Filtering. Least-Squares and Modeling: A Practical Handbook. Hoboken, New Jersey, John Wiley and Sons, 2011.

29. Semushin I.V., Tsyganova Yu.V., Kulikova M.V. et al. Adaptivnye sistemy fil'tracii, upravleniya i obnaruzheniya [Adaptive Systems of Filtering, Control and Betection]. Ulyanovsk, U1GU, 2011. (in Russian)

30. Semushin I.V. Vychislitel'nye metody algebry i otsenivaniya. Uchebnoe posobie [Computational Methods of Algebra and Estimation]. Ulyanovsk, U1STU Publishers, 2011. (in Russian)

31. Semushin I.V., Tsyganova Yu.V., Zakharov K.V. [Robust Filter Algorithms - Survey and New Results for Ship Navigation]. Informatsionnye tekhnologii i vychislitel'nye sistemy, 2013, no. 4, pp. 90-112. (in Russian)

32. Kulikova M.V., Taylor D.R. Stochastic Volatility Models for Exchange Rates and Their Estimation Using Quasi-Maximum-Likelihood Methods: an Application to the South African Rand. Journal of Applied Statistics, 2013, vol. 40, no. 3, pp. 495-507. DOI: 10.1080/02664763.2012.740791

33. Kulikova M.V., Tsyganova J.V., Semushin I.V. Adaptive Wave Filtering for Marine Vessels within UD-Based Algorithms. Proceedings of the ECC2016, European Control Conference, Aalborg, 2016, pp. 807-812. DOI: 10.1109/ECC.2016.7810388

34. Tsyganova Yu.V. [On the UD Filter Implementation Methods]. University Proceedings. Volga Region. Physical and Mathematical Sciences, 2013, no. 3, pp. 84-104. (in Russian)

35. Semushin I.V., Tsyganova Yu.V., Tsyganov A.V., Prokhorova E.F. Numerically Efficient UD Filter Based Channel Estimation for OFDM Wireless Communication Technology. Procedia Engineering, 2017, vol. 201, pp. 726-735. DOI: 10.1016/j.proeng.2017.09.597

36. Kulikova M.V., Tsyganova J.V. Improved Discrete-Time Kalman Filtering within Singular Value Decomposition. IET Control Theory and Applications, 2017, vol. 11, no. 5, pp. 2412-2418. DOI: 10.1049/iet-cta.2016.1282

37. Tsyganova Yu.V. Orthogonalized Block Methods for Parameter Identification of Discrete Linear Stochastic Systems: D.Sc. Thesis. Ulyanovsk, Ulyanovsk State University at Ulyanovsk, 2017.

38. Semushin I.V., Tsyganova J.V., Tsyganov A.V. Numerically Efficient LD-computations for the Auxiliary Performance Index Based Control Optimization under Uncertainties. Proceedings of the 17th IFAC Workshop on Control Applications of Optimization Yekaterinburg, Russia, October 15-19, 2018, Yekaterinburg, 2018, vol. 51, pp. 568-573.

39. D'Souza Ch., Zanetti R. Information Formulation of the UDU Kalman Filter. IEEE Transactions on Aerospace and Electronic Systems, 2018. DOI: 10.1109/TAES.2018.2850379

Received August 17, 2018