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

УДК [519.254+519.654]:629.5.05

Ю. В. Цыганова О МЕТОДАХ РЕАЛИЗАЦИИ Ш-ФИЛЬТРА1

Аннотация. Актуальность и цели. Фильтр Калмана является математическим инструментом, завоевавшим широкую популярность среди специалистов в области оценивания и управления. Но он имеет один существенный недостаток -неустойчивость по отношению к ошибкам машинного округления при его практической реализации на ЭВМ. Проблема ошибок машинного округления является неустранимой ввиду ограниченной разрядности представления вещественных чисел с плавающей запятой на ЭВМ. Однако можно существенно уменьшить влияние ошибок машинного округления в алгебраически эквивалентных реализациях фильтра Калмана, которые называют численно эффективными реализациями. Они основаны на различных математических методах факторизации ковариационных матриц ошибок оценок, участвующих в уравнениях фильтра. Целью работы является изучение основных методов построения UD-реализаций дискретного фильтра Калмана, обладающих улучшенными вычислительными свойствами по сравнению со стандартной реализацией фильтра Калмана, а также построение новой расширенной формы орто-гонализованного UD-фильтра, которая должна обладать следующими свойствами: устойчивость по отношению к ошибкам машинного округления; отсутствие операции извлечения квадратного корня; избавление от операции матричного обращения на каждой итерации алгоритма; компактность и удобство записи ортогонализованной формы UD-фильтра. Материалы и методы. Рассматриваются методы реализации UD-фильтров. Первой UD-реализацией фильтра Калмана является последовательный алгоритм Бирмана, а самыми современными являются ортогонализованные блочные алгоритмы. Подход к построению квадратно-корневых ортогонализованных блочных алгоритмов был предложен Кайлатом. В настоящей работе именно этот подход применяется для построения новой формы расширенного ортогонализованного UD-фильтра. Результаты. В работе изучены существующие к настоящему времени методы построения UD-фильтра. Наиболее эффективными в вычислительном плане и подходящими для реализации на современных вычислительных комплексах являются ортогонализованные формы UD-фильтра. Предложена новая форма расширенного ортогонализованного UD-фильтра, обладающая рядом преимуществ по сравнению с другими. Выводы. UD-алгоритмы в последовательной и матричной ортогонализованной формах являются эффективными в вычислительном плане реализациями дискретного фильтра Калмана. Их преимущества: 1) устойчивость по отношению к ошибкам машинного округления;

2) отсутствие операции извлечения квадратного корня; 3) избавление от операции матричного обращения на каждой итерации алгоритма; 4) компактность и удобство записи ортогонализованной формы UD-фильтра. Численные эксперименты показали работоспособность предложенной новой формы расширенного ортогонализованного UD-фильтра и ее устойчивость по отношению к ошибкам машинного округления на примерах плохо обусловленных задач.

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

1 Работа выполнена в рамках Государственных заданий Министерства образования и науки РФ (шифры 1.919.2011 и 6.3072.2011), поддержана Грантом РФФИ (шифр 13-01-97035).

Yu. V. Tsyganova ON THE UD FILTER IMPLEMENTATION METHODS

Abstract. Background. The Kalman filter (KF) is a mathematical tool that has won wide popularity among professionals in the field of estimation and control. However, it has a drawback - instability in respect to machine round off errors in its practical implementation on a computer. It is well known that the problem of machine round off errors is unavoidable due to the limited machine width of real floatingpoint numbers. However, one can significantly reduce the effect of round off errors in algebraically equivalent Kalman filter implementations (the so-called numerically efficient implementations). They are based on different mathematical factorization methods of the estimation error covariance matrices involved in the filter equations. The aim of the paper is to study the basic UD implementation methods of the discrete Kalman Filter with improved computational properties in comparison with the standard KF implementation, as well as the construction of the new extended orthogonal form of the UD filter which should have the following properties: robustness of computations against round off errors; lack of the square-root operation; deliverance from matrix inverse operation on each stage of the algorithm; compact and convenient orthogonal form of the UD filter. Materials and methods. The paper is based on the UD filter implementation methods. The first UD implementation of the KF was the Bierman’s sequential algorithm. Matrix orthogonal algorithms are the most modern numerically efficient KF implementations. The approach to the construction of the square-root orthogonal algorithms was proposed by Kailath and it is used in the paper to construct a new form of the extended orthogonal UD filter. Results. In the paper the existing implementation methods of UD filter have been investigated. The orthogonal forms of the UD filter are the most computationally efficient and suitable for implementation in modern computer systems. A new form of the extended orthogonal UD filter, which has a number of advantages in comparison with others, has been proposed. Conclusions. The sequential and orthogonal forms of the UD algorithms are the most numerically efficient implementations of the discrete Kalman filter. Their benefits are: 1) robustness of computations against round off errors; 2) lack of the square-root operation; 3) deliverance from matrix inverse operation on each stage of the algorithm; 4) compact and convenient orthogonal form of the UD filter. The numerical experiments have shown the efficiency of the proposed new form of extended orthogonal UD filter and its computations robustness against round off errors while dealing with ill-conditioned tasks.

Key words: stochastic discrete linear systems, optimal discrete filtering, the Kalman filter, robust implementations of the discrete filter, the UD filter.

Введение

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

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

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

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

В дальнейшем эти и другие подобные методы стали отправной точкой для развития целого направления в теории линейного оценивания и управления, которое продолжает активно развиваться и в настоящее время. Существует огромное количество научных публикаций по данной теме. Здесь отметим только хорошо известные автору монографии, содержащие описание численно эффективных реализаций ФК: зарубежные [3-6] и отечественные [7-9].

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

UpDpUp , где U - верхняя треугольная матрица с единицами на диагонали,

D - диагональная матрица (для любой квадратной положительно определенной матрицы такое представление можно получить в результате модифицированного разложения Холецкого [3, 6]). Первой UD-реализацией ФК является последовательный алгоритм Дж. Бирмана [3]. Бирман не только доказал вычислительную эффективность UD-фильтра, но и показал, что при соответствующей программной реализации его алгоритм не сложнее, чем стандартный алгоритм Калмана [10].

Из всех численно эффективных реализаций ФК самыми современными являются ортогонализованные блочные алгоритмы. Подход к построению квадратно-корневых ортогонализованных блочных алгоритмов был предложен Т. Кайлатом (подробное описание см. в [5]). В настоящей работе именно этот подход применяется для построения новой формы расширенного орто-гонализованного UD-фильтра.

1. Постановка задачи

Рассмотрим стохастическую дискретную линейную систему:

4+1 = фкХк + Gkwk, wk е N(0,^X (1)

ч = Hkxk + ч, че N (0, Ч X k ^ 0 (2)

где Хк е RП - вектор состояния системы, Zk е Rт - доступный вектор измерений. Последовательности ^0,Wl,...} и {^,У2,...} - независимые нормально распределенные белые последовательности шумов с нулевыми

средними и ковариационными матрицами Qk ^ 0 и Я >0 соответственно,

є RпХі. Последовательности } и [у.} не зависят от начального состояния Хо є N(хо,По). Другими словами, выполняется следующее условие:

wk 'Qk 0 " S 0 Ski kk 0 0 '

E vk [Wj vT xl 1] = 0 Rk _ 0

x0 o 0 П0 x0 _

где - символ Кронекера; E[•] - символ математического ожидания.

Для решения задач рекуррентного оценивания и адаптивной фильтрации [5, 7, 11], параметрической идентификации [12], стохастического управления [4, 13, 14] и многих других применяется широко известный со второй половины прошлого века математический инструмент - дискретный фильтр Калмана [1].

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

Алгоритм 1 (двухстадийная форма стандартной реализации

дискретного ФК).

0. НАЧАЛЬНЫЕ ДАННЫЕ: х0|0 = х0 и _Р0|0 = По , П0 > 0.

1. РЕКУРРЕНТНО ОБНОВЛЯТЬ ВЕЛИЧИНЫ (к > 0 ).

1.1. НА ЭТАПЕ ФИЛЬТРАЦИИ:

Кк = ркнк К-к, Кв,к = Кк + нкркнк (коэффициент КалманаХ хк|к = Хк + Ккек, ек = 2к -НкХк (°ценкаХ

Рщ = Рк - КкНкРк (матрица ковариации ошибки).

1.2. НА ЭТАПЕ ЭКСТРАПОЛЯЦИИ:

Хк+1 = фкхк |к (оценкаХ

Рк+1 = ФкРк|кФк +^0^ (матрица ковариации ошибки), (3)

где %к = -Хк|к-1 - линейная среднеквадратичная экстраполированная оценка вектора состояния Хк, вычисленная по доступным измерениям { ^,..., Zk — } ,

Рк = Рк\к-1 = Щ(Хк - Хк\к-1)(Хк - Хк\к-1)Т] -

ковариационная матрица ошибки оценивания на этапе экстраполяции; Хк|к -линейная среднеквадратичная отфильтрованная по текущему измерению Zk

T

оценка вектора состояния, Pk|k = E[(xk - xk|k)(xk - xk|k) ] - ковариационная

матрица ошибки оценивания на этапе фильтрации.

Если нас интересует только экстраполированная оценка вектора состояния и соответствующая ей ковариационная матрица ошибки оценивания, то алгоритм 1 можно переписать в удобной одностадийной форме [5].

Алгоритм 2 (одностадийная форма стандартной реализации

дискретного ФК).

0. НАЧАЛЬНЫЕ ДАННЫЕ: Jcojo = xo и Po|o = П0 , П0 >0 .

1. РЕКУРРЕНТНО ОБНОВЛЯТЬ ВЕЛИЧИНЫ (k > 0):

Kpk = ^PkHk Re-k, Rek =Rk + HkPkHk (матрица обратной связи), (4)

Pk+1 = фkpkфT + GkQkGk - Kp,kRe,kKP,k (матрица ковариации °шибкиХ (5)

xk+1 = ФЛ + Kp,kek, ek = zk -Hkxk (оценка). (б)

Далее в работе будем использовать следующие обозначения: для

-1 T T -1 -T

невырожденной квадратной матрицы (A ) = (A ) = A ; через D

обозначим диагональную матрицу, через U и L - соответственно верхнюю и нижнюю треугольную матрицу с единицами на диагонали. Используем модифицированное UD-разложение Холецкого [З, б] для положительно

определенной квадратной матрицы P = UpDpUp .

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

- устойчивость по отношению к ошибкам машинного округления;

- отсутствие операции извлечения квадратного корня;

- избавление от операции матричного обращения на каждой итерации алгоритма;

- компактность и удобство записи ортогонализованной формы UD-фильтра.

2. Методы реализации UD-фильтра

Перечислим базовые идеи, на которых построены различные формы UD-реализации дискретного ФК:

1. Дж. Бирман (G. Bierman) [З] предложил использовать UD-факторы ковариационной матрицы ошибки оценивания вектора состояния для обновления разностного уравнения Риккати на этапе обработки измерений в фильтре Калмана.

2. К. Торнтон (C. Thornton) [15] применила алгоритм взвешенной орто-гонализации Грама - Шмидта для вычисления UD-факторов ковариационной матрицы на этапе экстраполяции в фильтре Калмана.

3. Т. Кайлат (T. Kailath) [16] построил новый класс алгоритмов для реализации дискретного ФК - квадратно-корневые ортогонализованные мат-

ричные алгоритмы (array algorithms). Данные алгоритмы базируются на применении ортогональных преобразований к блочной матрице, называемой pre-array. В результате получается новая блочная матрица, называемая postarray, с готовыми результатами вычислений. Отличительной чертой в таких алгоритмах является удобная форма представления данных для реализации на ЭВМ, ориентированность на параллельные вычисления [17, 18]. Кроме того, ортогональные преобразования обладают численной устойчивостью по отношению к ошибкам машинного округления.

4. Т. Кайлат [17] предложил ортогонализованный алгоритм этапа обработки измерения в ФК, основанный на быстрых вращениях Гивенса. Было показано, что данный алгоритм эквивалентен алгоритму Бирмана для случая скалярных измерений.

5. В книге [б] авторами сформулирован одностадийный ортогонализо-ванный алгоритм для UD-фильтра, основанный на модифицированной взвешенной ортогонализации Грама - Шмидта.

3.1. Последовательный алгоритм Бирмана

Пусть на этапе фильтрации в алгоритме 1 используются разложения

Pk = Up Dp ul и Pkk = Up Dp ul . k pk pk pk k|k Pkk pk\k pk\k

Предположим, что элементы матрицы Dp являются положительными,

k

а матрица Rk - диагональная, т.е. Rk =Diag[rk(1),...,rk(m)]. Тогда этап фильтрации алгоритма 1 эквивалентен следующему алгоритму:

Алгоритм 3 (последовательный алгоритм Бирмана).

0. НАЧАЛЬНОЕ ПРИСВАИВАНИЕ: U = Up , D = DP , X = Xk.

kk

1. m -КРАТНОЕ ПОВТОРЕНИЕ ПРОЦЕДУРЫ СКАЛЯРНОГО ОБНОВЛЕНИЯ.

Для j = 1,2, ..,m выполнять:

1.1. Вычислить векторы

f = fi,/2,., fn ]T= UTh;

v = fv2,.,vn]T = Df .

1.2. Задать начальные значения а' = r + vi fi, di = dir / аі;

K = [vi |0.0]T .

1.3. Для i = 2,..,n выполнять:

начало

а := а'+ vif; Y := 1/ а;

di := dj-а' y; X := -fiY;

Hi := Hij + ^K; K := K+vHi;

а' := а.

конец

1.4. Вычислить векторы V := у(г - И х); х:= х + КV с экстраполяцией между повторениями: ) = и; Г) = Г ; х = х.

2. ЗАВЕРШАЮЩЕЕ ПРИСВАИВАНИЕ:

up = U: pk\k

Dp = D;

pk|k

T

хк|к - х •

Здесь к - у-й столбец матрицы Нк ; г - у-й элемент вектора ^ ; г - у-й элемент Гу (к) диагональной матрицы ковариаций шума измерений Як, у = 1,..,т - номер скалярного измерения в составе вектора измерений ^ в момент времени к • Доказательство алгоритма можно найти в [3, с. 77-80].

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

3.2. Алгоритм Торнтон на базе MWGS-ортогонализации

Предположим, что на этапе экстраполяции доступны иБ-факторы

гк|к~-к|к~рк|к “ -^-^ак

Т Т

ковариационных матриц Рщ = иРщЕрщир и Як = ■ Тогда

уравнение (3) принимает вид

Pk+1 = ФkUp., Dp„ U

JkuP,

T

k|k Pk\k pk\k

ФT + GkJ

kUgkDgkuQkGT

Представим матрицу Рк+1 как

Pk+1 =[Ф kUpk\k \GkJgk] DiagiDpk\k ,Dgk

rk\k

rk\k'

u

T

p.

k\k

T

gk

ФT

ux GT

т.е. Pk+1 = WkDkWi

T

k^k^k

Wk =[Фk|kUPk\k \GkUgk]:

Dk =Diag{Dp , Dgk } = Diag {di(k),...,dM (k)},

wiT(k)

wl(k)

где Щ - матрица размера п X (п + 5) и М = п + 5 , 5 - размер вектора шума Wfc в уравнении состояния (1). Факторы )р^ ^ и 1 вычисляются

k+1

методом модифицированной взвешенной ортогонализации Грама - Шмидта ([3], теорема VI.4.1) на этапе экстраполяции (описание метода MWGS приведено в приложении). Сформулируем алгоритм Торнтон.

Алгоритм 4 (алгоритм Торнтон для этапа экстраполяции в ФК).

0. НАЧАЛЬНЫЕ ДАННЫЕ: ^, {Ц^, Ор^}, {и^, }.

1. ВЫЧИСЛИТЬ ФАКТОРЫ {иР , БР }:

1 Рк+1 Рк+г

1.1. Сформировать блочные матрицы

°к =°іа§{БРк|к В2к } Щ = [Фк|киРк,к \GkUQk].

1.2. Для пары {Щ,Вк} провести MWGS-ортогонализацию строк матрицы Щ относительно весовой матрицы Бк и получить в результате тру {иРк+г Орк+1}.

1.3. Вычислить оценку вектора состояния: Хк+1 = ФкХк|к .

Доказательство алгоритма Торнтон можно найти в [3, с. 127-128]. Этот алгоритм также не содержит операции извлечения квадратного корня, а однородная структура алгоритма подходит для реализации на современных вычислительных комплексах.

3.3. Одностадийная ортогонализованная форма иБ-фильтра

Т. Кайлат [17] предложил ортогонализованный алгоритм этапа обработки измерения в ФК, основанный на быстрых вращениях Гивенса [19]. Алгоритм был сформулирован для ЬБ-факторов ковариационных матриц.

Предположим, что на этапе обработки измерения доступны ЬБ-

и рк = . Тогда

факторы ковариационных матриц Pk = Lp Dp L

kpk

LD-факторы Lpkk и Dp^k ковариационной матрицы Pk|k можно получить

k|k

в результате матричного ортогонального преобразования

1 L k HkLpk - i D k О " = Q /2 Г lr о 1 Re,k \dr Re,k 0 '

О I L О Й* D _Kk\k Lpk\k _ 0 I D

1/2

, (7)

где 2 - матрица ортогонального преобразования; Ьр ^ и Бр к

- ЬБ-факторы

ковариационной матрицы невязки измерений Рек; Кк - матрица Калмана.

Если вместо 2 использовать Яв матрицу преобразования на основе быстрых вращений Гивенса, не содержащих операцию извлечения квадратного корня, то соотношение (7) можно переписать в виде

1 L k HkLpk - qd = Г lr о 1 Re,k

0 I Lpk L _Kk\k Lpk\k I

В работе [17] было показано, что данный алгоритм эквивалентен ЬБ-реализации алгоритма Бирмана для случая скалярных измерений.

На базе алгоритмов Т. Кайлата в книге [6, с. 262] сформулирован одностадийный ортогонализованный алгоритм для иБ-фильтра. Вместо быстрых вращений Гивенса применяется модифицированная взвешенная

ортогонализация Грама - Шмидта. Запишем указанный алгоритм для системы (1), (2).

Алгоритм 5 (одностадийная ортогонализованная форма ИБ-фильтра).

0. НАЧАЛЬНЫЕ ДАННЫЕ: х0 =0 и Р0 = П0 > 0 . Вычислить ИБ-

факт°ры Холецкого: {ип0, ^п0 }, {UR0, °я0 } , {ие0, ^е0 } .

1. РЕКУРРЕНТНО ОБНОВЛЯТЬ {иР ,ВР } и хк+1 (к > 0):

£+1 к+1

1.1. Для {ирк, Брк }, {иКк, 0Кк }, {и^к, Бдк } сформировать блочные матрицы

Dk =Dlag{DQk , DPk , DRk }, Ak

GkUQk ФkUPk 0

0 HkUPk URk

1.2. Для пары ^к, Dk} провести MWGS-ортогонализацию столбцов матрицы Ak относительно весовой матрицы Dk и получить в результате пару блочных матриц ^к, Djk} :

Dk =Dag{DPk+1, DRek}, Ak =

UP Kp kUR

Pk+l p,k Rek

О UR

Re,k

таких, что

А-І “ AкТвк и AkТDkAk “ Aк^ккТ)1,

Б т>(п+т+і)х(п+т+і) д _ лз(п+т+і)х(п+т) п _ о(п+т+в)х(п+т)

к є R , A к є R , вк є R -

матрица MWGS-преобразования, приводящего к верхней треугольной блочной матрице AкТ є R(п+т)х(п+т) и диагональной матрице DkТ є R(п+т)х(п+т)

1.3. Вычислить оценку вектора состояния:

Хк+1 = фкХк + (Кр,киР к ) ек, где ек = и-\ ек, ек = Ч - Нхк. (8)

е,к ре,к

Отметим, что величина Кр ^р ^ (в круглых скобках в (8)) получена

непосредственно из AкТ .

Алгоритм 5 позволяет вычислить те же величины, что и алгоритм 2, но, в отличие от последнего, является устойчивым по отношению к ошибкам машинного округления. Вместе с тем он обладает следующими недостатками:

1) вычисление оценки вектора состояния по уравнению (8) происходит отдельно от обработки остальных данных, которые входят в блочные матрицы и обрабатываются с помощью MWGS-алгоритма. Таким образом, алгоритм имеет неоднородную структуру;

2) в уравнении (8) присутствует операция матричного обращения, которая может ухудшить качество оценок вектора состояния.

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

4. Расширенная ортогонализованная форма UD-фильтра

Идея расширенного ортогонализованного алгоритма заключается в том, что все величины дискретного фильтра помещаются в блочную матрицу и обрабатываются с помощью ортогональных преобразований. Впервые эта идея была реализована в [20] для построения квадратно-корневого информационного фильтра (SRIF - square-root information filter). Для ковариационной формы ФК расширенные квадратно-корневые ортогонализованные алгоритмы построены в [16].

Рассмотрим снова представление ковариационной матрицы в виде

pk=UPkDPkUT.

Введем обозначения:

4 - (UPkDPk ) 1 Хк ,

Ьк = -(иЯеЛ,к У1вк ■

Сформулируем алгоритм, реализующий расширенную ортогонализованную форму ИБ-фильтра.

Алгоритм 6 (расширенная ортогонализованная форма ИБ-фильтра).

0. НАЧАЛЬНЫЕ ДАННЫЕ: £0 - 0 и Р0 = П0 > 0 . Вычислить ИБ-

факт°ры Холецкого: {ип0, £Л0}, (UR0, °я0} , {UQ0, °(20 } ■

1. РЕКУРРЕНТНО ОБНОВЛЯТЬ (и^+1,Орк+1} и гк+1 (к >0):

1.1. Для (ирк, Врк }, (иКк, DRk }, {^к, }, 2к сформировать

блочные матрицы

% = Б*а§ (DQk, DPk, ^к },

0 T zk -zT(URkDRk )-T

II Tk A GkUQk ф^ 0

0 HkUPk URk

(9)

1.2. Для пары {A^,D^} провести MWGS-ортогонализацию столбцов матрицы Ak относительно весовой матрицы Dk и получить в результате

А А

пару блочных матриц ^^, Dk } :

dI — Diagfl, DPk+1-dr },

(*)

0 Upk

о

~T zk+1

k+1 о

bk

KpkU\k

U

R

e,k

(1о)

таких, что

и

Ab Ak1 BkT

A k T Dk A k — A k IDk I(A k I)T,

D d( n+m+s)x(n+m+s) » D( n+m+s )x(1+n+m) R D(n+m+s)x(1+n+m)

k є R - Ak є R - Bk є R —

матрица MWGS-преобразования, приводящего к верхней треугольной

и диагональной матрице

1; (*) - величина, не представляющая для нас интерес. Замечание. В любой момент времени к доступна оценка вектора состояния

(11)

блочной матрице Ak1 є R(1+n+m)x(1+n+m)

Dk1 є R(1+n+m)x(1+n+m) • (*)

xk —(UPkDPk) zk.

Т

Доказательство алгоритма 6. Представим произведение Aк DkA

в виде

(*) A2T1 A31

A21 A22 A,2

A31 A32 A33

Ak DkAk —

Определим элементы блочной матрицы (12):

А21 = ф киРкВркїк = ф%;

А22 = °кидк»дкидкск +ФкиркОркиТк Ф^ = + ФкРк Ф^

А31 = нкиРквРк^к -Ч = нкХк -Ч = -ек;

(12)

lkUpkUpkUPk

T

T

A32 — HkUp Dp U* Фk — HkPkфk;

A33 — HkUp Dp UT + Ur, Dr, UT — HkPkHj + Rk.

IhUp Up U n + Up Up U n —Пь-Гь-Пь-

k Pk Pk Pk Rk Rk Rk k k k

А А А Т

Аналогично представим произведение Л к Dk (А к) в виде

A1 Dlk (A! )T —

(*) A2T1 A,Tl

A21 A22 A32

A31 A32 A33

(13)

Определим элементы блочной матрицы (13):

а21 = ^+1 Ърк+14+1 + крЛиКекоКекьк =

= -КркиК 0 (ЧА,к)-1 ек = ^ -Кркек;

A22 = Up DP uT + KpkUR Dr uT Klk =

22 Pk+1 Pk+1 Pk+1 p’k Re,k Re,k Rek P,k

T

Pk+1 + Kp,kRe,kKp,k;

A31= URekDRekbk = -ek;

A32= = Re,kK.Tk;

A33 - UR kDR kUi^ - Re,k•

e,k e,k Re,k

Сравнивая левую (12) и правую (13) части равенства т ^ ^ ^ т

А к Бк А к = А к Бк (А к) , приходим к соотношениям:

А21 _ А21 ^ Фхк = Хк+1 — Кр,кек ^ (6);

А22 = А22 ^ СкЯкбк + ФкРк Фт = Рк+1 +Кр,кРе,кКр,к ^ (5);

А31 = А31; А32 = А32;

А33 = А33 ^ НкРкНк + Рк = Ре,к ^ (4)

что и доказывает справедливость алгоритма 6.

Построенный алгоритм 6 обладает заявленными в разделе 2 свойствами. Он имеет компактную и удобную форму записи, позволяющую единообразно обрабатывать все данные. Такая форма иБ-фильтра позволяет предельно просто реализовать его в виде программного кода на ЭВМ. В алгоритме отсутствует операция извлечения квадратного корня, а вычисление оценки вектора состояния (11) теперь не требует операции матричного обращения.

Замечание 2. Для того чтобы избавиться от матричного обращения при

т —т

вычислении элемента блочной матрицы — Zk (иркБрк) в (10), можно

* т *

заменить его на —(Zk) , где Zk - декоррелированные измерения (см. [9, с. 240]).

5. Численные примеры

Рассмотрим прикладную задачу [9] получения оценок высоты и вертикальной скорости летательного аппарата (ЛА) по показаниям двух приборов: инерциального датчика вертикального ускорения ау (^) и барометрического

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

Введем обозначения физических величин, участвующих в математической модели движения объекта (ЛА) по высоте и в модели наблюдения за этим

движением. В соответствии с принципиальной моделью - вторым законом Ньютона, движение центра масс объекта характеризуем тремя величинами: У = у(() - вертикальная координата (высота над Землей), у у = Уу (() -

вертикальная скорость, ау = ау (() - вертикальное ускорение.

Будем рассматривать случай движения с постоянной силой тяги двигателя ЛА. В этом случае следует считать, что ускорение ЛА тоже постоянно (пренебрегая изменением массы ЛА), однако оно неизвестно, т.е. подлежит оцениванию по показаниям приборов. В роли оценивателя наряду со стандартным фильтром Калмана используем рассмотренные выше ИБ-реализации.

Полученные показания приборов (измеренные величины) обозначим: za = za (() - измеренное ускорение а((), Zh = Zh (() - измеренная

барометрическая высота й(().

Считывание этих показаний происходит в дискретные моменты времени с аддитивными погрешностями у^-) и У2 ((-) соответственно. Это означает, что za (() = а((-) + У1 ((-), Zh ((-) = И((-) + У2((-). Погрешности датчиков между собой независимы и в отдельные моменты времени являются

гауссовыми случайными величинами с нулевыми средними значениями и

22

с постоянными дисперсиями ©1 , ©2 соответственно. Истинная высота у((), на которой находится ЛА, воздействует на барометрический датчик (БД) высоты так, что барометрическая высота й(() подчиняется уравнению

т—И(() + И(( ) = у(().

Запишем уравнения состояния обобщенного динамического объекта, включающего ЛА и БД. Обозначим вектор состояния как

х = х(() = [Х1 = у(()|Х2 г Уу (()|хз = ау(()|Х4 = И(г)]т

и вектор измерения как

z = z (() =[ Z1 = za(()| z2 = Ч(()]т т

с погрешностью у(() = У(()|у2 (()] . Из принятых предположений получаем непрерывную модель состояния

—x(t) — Fx(t), F — dt

и дискретную модель наблюдения

Ч = Нхк + Ч, н =

0 1 0 0

0 0 1 0

0 0 0 0

a 0 0 —а

1

а--, т

0 0 10 0 0 0 1

Далее, переходя к дискретной модели состояния, запишем переходную матрицу состояния на интервале времени (:

—at

1

0

0

t2/2 t 1

al(t) a2(t) a,(t)

0

0

0

a4(t)

al(t) — 1 — e a2(t) — (at — 1 + e—at)/a,

a,(t) — (1 — at + a1(t) — e—at.

(at )2

-e—at)/a2,

Зададим т* - постоянный интервал выборки, т.е. темп считывания данных с датчиков; при этом считаем, что у = т* /т<1, например, у = 1/10 или меньше. Подстановкой ( = т* в Ф(() определим постоянную Ф = Ф(т*) -переходную матрицу состояния дискретной модели

Г Хк+1= фхк + вмк, ^ = Нхк + ук ]

Модель ЛA + БД:

G — [0 1 0 0] — const

где Gwk - «шумовая» составляющая с независимым дискретным белым шумом Wfc е N (0, Q = qT s = const), q - коэффициент диффузии соответствующего винеровского процесса, м2 / с5 (учитываем нестационарность силы тяги двигателя ЛА).

Для заданной модели «ЛА+БД» проведем вычислительные эксперименты. Экспериментальные условия приведены в табл. 1.

Таблица 1

Экспериментальные условия для модели «ЛЛ+БД»

№ Параметр Вариант исходных данных

1 2 3 4 5 б

1 Диффузия д, м2/с5 3000 340 400 300 3500 3350

2 Постоянная времени т, с 0,05 0,б5 0,80 0,90 0,10 0,12

3 Дисперсия о2, (м/с2)2 1,00 2,25 4,00 б,25 б,00 5,50

4 Дисперсия ©2 , м 40 30 25 35 45 50

5 Дисперсия Рп (0), м2 10 20 30 40 50 б0

б Дисперсия Р22 (0), (м/с)2 б0 50 40 30 20 10

7 Дисперсия Р33 (0), (м/с2)2 15 20 25 30 35 40

В Дисперсия Р44 (0), м2 45 40 35 25 30 15

Выберем количество измерений N равным 100. Сначала экспериментально оценим разницу в вычислениях по следующим алгоритмам:

1) 8КБ - стандартная реализация дискретного ФК (алгоритм 2);

2) ТВБ - дискретный фильтр Торнтон - Бирмана (алгоритмы 3 и 4);

3) еИБ - расширенный ортогонализованный ИБ-фильтр (алгоритм 6). Обозначим через Аа_в (X) разницу в вычислениях оценки вектора

состояния по алгоритмам А и В; при этом

^АВ (X) = max II Xk - xk ||то,

1<k < N

где Xk и Xk - оценки вектора состояния Xk , вычисленные по алгоритмам А и В в равных экспериментальных условиях. Аналогично обозначим через А А,В (P) разницу в вычислениях ковариационной матрицы ошибки

оценивания по алгоритмам А и В; при этом

ДА,В(P)= max 11 Pk -Pk !L,

1<k < N

А b

где Pk и Pk - ковариационные матрицы ошибки оценивания, вычисленные по алгоритмам А и В в равных экспериментальных условиях. Результаты вычислительных экспериментов приведены в табл. 2 и 3. Анализируя полученные данные, можно сделать вывод о том, что разница в вычислениях по алгоритмам SKF, TBF и eUD не превышает 10-12, т.е. в условиях рассмотренной прикладной задачи данные алгоритмы практически эквивалентны. Этот факт подтверждает работоспособность и правильность построения алгоритма eUD.

Таблица 2

Разница в вычислениях оценки вектора состояния x для модели «ЛА+БД»

Алгоритм А A,B (X)

A B 1 2 3 4 5 б

SKF TBF 2,13 • 10-14 б,39 • 10-14 1,99 • 10-13 5,8б • 10-14 3,55 • 10-14 7,11 • 10-14

SKF eUD 2,13 • 10-14 3,б2 • 10-13 7,18 • 10-13 2,42 • 10-13 2,27 • 10-13 1,28 • 10-13

eUD TBF 1,73 • 10-14 3,87 • 10-13 7,39 • 10-13 2,49 • 10-13 2,20 • 10-13 1,71 • 10-13

Таблица 3

Разница в вычислениях ковариационной матрицы ошибки оценивания Р для модели «ЛА+БД»

Алгоритм А A,B (P)

A B 1 2 3 4 5 б

SKF TBF 1,02 • 10-12 7,39 • 10-13 2,27 • 10-13 3,98 • 10-13 2,27 • 10-13 2,05 • 10-12

SKF eUD 9,09 • 10-13 7,39 • 10-13 2,27 • 10-13 3,98 • 10-13 5,б8 • 10-13 2,05 • 10-12

eUD TBF 7,9б • 10-13 2,27 • 10-13 3,41 • 10-13 2,27 • 10-13 б,82 • 10-13 1,14 • 10-12

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

Рассмотрим задачу оптимальной дискретной фильтрации. Матрицы-параметры системы (1), (2) заданы следующим образом:

Ф k = /з, Qk =0, Gk =0, Rk = б2 /2, Hk =

11 1 1 1 1 + б

2

начальные условия х0 є N(0, /3), где 5 < £окр , но 5 > £окр ; £окр - параметр

окр

машинного округления, т.е. єокр +1Ф1, но єокр /2 +1 = 1. Требуется вычислить значение матрицы ковариации ошибки оценивания Р после обработки одного измерения, т.е. Р[|1. Данный тестовый пример показывает, как влияют

ошибки машинного округления на работоспособность ковариационных фильтров.

Легко видеть, что ранг матрицы Нк равен 2. Можно показать, что при

малых значениях 5 матрица ковариации невязки измерений Яе 1 является

плохо обусловленной, а при 5 ^ £окр в результате округления она имеет

ранг, равный 1. В работе [21] проведен строгий теоретический анализ влияния ошибок округления на некоторые из реализаций дискретного фильтра Калмана. Авторами доказано, что число обусловленности матрицы ковариации невязки измерений Яе к является определяющим фактором,

влияющим на численное поведение фильтров ковариационного типа.

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

ТББ, еИБ и квадратно-корневому алгоритму Поттера [3, с. 22-24; 7, с. 79]. Для каждого значения 5 будем вычислять максимальную относительную ошибку среди элементов матрицы Рщ (обозначим АР1), сравнивая точное

решение задачи с численным, полученным с использованием ЭВМ.

Результаты экспериментов представлены в виде графиков на рис. 1. Для каждой реализации дискретного фильтра изображена зависимость АР от параметра обусловленности 5, где 5 ^ £окр .

Численные расчеты проводились на персональном компьютере, обеспечивающем представление чисел с относительной погрешностью 10-16. Коды программ написаны на языке МАТЬАБ, где параметр машинного округления £окр хранится в переменной ер8 со значением £окр = ер8/2 =

= 1,1102•10_16 . Проведем расчеты для 5є [10-9 ер82/3,109ер82/3].

Проанализируем результаты эксперимента. Легко видеть, что стандартный алгоритм Калмана (сплошная линия с маркером *) для рассмотренной задачи при 5 ^ £окр быстро теряет точность вычислений. Так, например, при _______8

5 = 10 стандартная реализация фильтра обеспечивает вычисление Рщ

с точностью 10-1, в то время как ИБ-реализации (сплошные линии

с маркерами X и о) обеспечивают точность 10 9. Таким образом, уже при

5 = 10 8 и меньше алгоритм 8КБ не обеспечивает даже одной правильной цифры после запятой, т.е. теряет работоспособность. ИБ-реализациям удается поддерживать приемлемую точность вычислений вплоть до граничного значения 5 = £окр .

5

Рис. 1. Сравнение точности вычислений по четырем алгоритмам: стандартный алгоритм Калмана (SKF), UD-реализация Торнтон - Бирмана (TBF), расширенная UD-реализация (eUD), квадратно-корневой алгоритм Поттера (Potter)

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

Заключение

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

Основная цель дальнейших исследований заключается в построении адаптивного расширенного ортогонализованного UD-фильтра для эффективного в вычислительном плане решения задачи параметрической идентификации дискретных моделей линейных стохастических систем. Подобные методы построения адаптивных ортогонализованных фильтров разработаны в [22-24].

Приложение

Модифицированная взвешенная ортогонализация Грама - Шмидта Алгоритм MWGS-ортогонализации был построен А. Бьерк (А. Bjбrck [25]). Доказано [26], что данный алгоритм эффективнее в вычислительном плане, чем классическая ортогонализация Грама - Шмидта, и точность вычислений по нему сравнима с точностью вычислений по известному алгоритму триангуляризации Хаусхолдера.

В алгоритме взвешенной ортогонализации рассматривается ортогональность т -векторов Ь и Ъу относительно весовой матрицы О^, :

ЪТ О Ъ _ ГР/ >0, если 1 _ ■,

Ъ1 °М>Ъ] | о ■ Ф ■

[ 0, если / Ф у.

Процесс MWGS-ортогонализации набора п линейно независимых т -векторов 0{, а2, ..., ат относительно весовой матрицы О^, приводит

Т Т

к верхней треугольной (п X п)-матрице и такой, что А _ ив , т.е.

Г t 1 a1 "1 u12 •• u1n

T a2 = 0 1 . u2n

T _ an _ 0 0 .. 1

bi

b2T

где вектора bi есть столбцы матрицы MWGS-преобразования B и

B DWB = Diag{pf } = Dp.

1<i<n

Список литературы

1. Kalman, R. E. A New Approach to Linear Filtering and Prediction Problems / R. E. Kalman. // Trans. of the ASME-Journal of Basic Engineering. - 1960. - Vol. 82, Series D. - P. 35-45.

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

3. Bierman, G. J. Factorization methods for discrete sequential estimation / G. J. Bierman : New York : Academic Press, 1977. - 255 p.

4. Maybeck, P. S. Stochastic models, estimation and control, Vol. 1. / P. S. Maybeck. -New York : Academic Press, Inc., 1979. - 423 p.

5. Kailath, T. Linear estimation / T. Kailath, A. H. Sayed, B. Hassibi. - New Jersey : Prentice Hall, 2000. - 854 p.

6. Grewal, M. S. Kalman filtering: theory and practice / M. S. Grewal, A. P. Andrews. -New Jersey : Prentice-Hall, 2001. - 410 p.

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

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

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

10. Bierman, G. J. Numerical Comparison of Kalman Filter Algorithms : Orbit Determination Case Study / G. J. Bierman, C. L. Thornton // Automatica. - 1977. -Vol. 13. - P. 23-35.

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

12. Льюнг, Л. Идентификация систем. Теория для пользователя : пер. с англ. / Л. Льюнг ; под. ред. Я. З. Цыпкина. - М. : Наука, 1991. - 432 c.

13. Саридис, Дж. Самоорганизующиеся стохастические системы управления : пер. с англ. / Дж. Саридис ; под. ред. Я. З. Цыпкина. - М. : Наука, 1980. - 400 с.

14. Mosca, E. Optimal, predictive and adaptive control / E. Mosca. - New Jersey : Prentice-Hall, Inc., 1995. - 477 p.

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

16. Park, P. New square-root algorithms for Kalman filtering / P. Park, T. Kailath // IEEE Trans. Automat. Control. - 1995. - V. 40, № 5. - P. 895-899.

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

18. Hotop, H.-J. New Kalman filter algorithms based on orthogonal transformations for serial and vector computers / H.-J. Hotop // Parallel Computing. - 1989. - № 12. -P. 233-247.

19. Gentleman, W. M. Least squares computations by Givens transformations without square roots / W. M. Gentleman // J. Inst. Math. Appl. - 1973. - № 12. - P. 329-336.

20. Dyer, P. Extension of square-root filtering to include process noise / P. Dyer,

S. McReynolds // J. Optim. Theory Appl. - 1969. - № 3. - P. 444-459.

21. Verhaegen, M. Numerical aspects of different Kalman filter implementations / M. Verhaegen, P. Van Dooren // IEEE Trans. Automat. Contr. - 1986. - Vol. 31, № 10. -P. 907-917.

22. Цыганова, Ю. В. Вычисление градиента вспомогательного функционала качества в задаче параметрической идентификации стохастических систем / Ю. В. Цыганова // Автоматика и телемеханика. - 2011. - № 9. - C. 142-160.

23. Цыганова, Ю. В. Об эффективных методах параметрической идентификации линейных дискретных стохастических систем / Ю. В. Цыганова, М. В. Куликова // Автоматика и телемеханика. - 2012. - № 6. - C. 34-51.

24. Tsyganova, J. V. State sensitivity evaluation within UD based array covariance filters / J. V. Tsyganova, M. V. Kulikova // IEEE Transactions on Automatic Control. -2013. - Vol. 58, № 11. - Р. 2938-2944.

25. Bjorck, A. Solving least squares problems by orthogonalization / A. Bjorck // BIT. -1967. - Vol. 7. - P. 1-21.

26. Jordan, T. L. Experiments on error growth associated with some linear least-squares procedures / T. L. Jordan // Math. Comp. - 1968. - Vol. 22. - P. 579-588.

References

1. Kalman R. E. Trans. of the ASME-Journal of Basic Engineering. 1960, vol. 82, Series D, pp. 35-45.

2. Potter J. E., Stern R. G. Proceedings of 1963 AIAA Guidance and Control Conference. New York: AIAA, 1963.

3. Bierman G. J. Factorization methods for discrete sequential estimation. New York: Academic Press, 1977, 255 p.

4. Maybeck P. S. Stochastic models, estimation and control, Vol. 1. New York: Academic Press, Inc., 1979, 423 p.

5. Kailath T., Sayed A. H., Hassibi B. Linear estimation. New Jersey: Prentice Hall, 2000, 854 p.

6. Grewal M. S., Andrews A. P. Kalman filtering: theory and practice. New Jersey: Prentice-Hall, 2001, 410 p.

7. Ogarkov M. A. Metody statisticheskogo otsenivaniya parametrov sluchaynykh protsessov [Methods of statistical estimation of stochastic process parameters]. Moscow: Energoatomizdat, 1990, 208 p.

8. Semushin I. V., Tsyganova Yu. V., Kulikova M. V. et al. Adaptivnye sistemy fil’tratsii, upravleniya i obnaruzheniya: kollektivnaya monografiya [Adaptive systems of filtration, control and detection: collaborative monograph]. Ulyanovsk: Izd-vo UlGU, 2011, 298 p.

9. Semushin I. V. Vychislitel’nye metody algebry i otsenivaniya: uchebnoe posobie [Computational methods of algebra and estimation: tutorial]. Ulyanovsk: UlGTU, 2011, 366 p.

10. Bierman G. J., Thornton C. L. Automatica. 1977, vol. 13, pp. 23-35.

11. Fomin V. N. Rekurrentnoe otsenivanie i adaptivnaya fil’tratsiya [Recurrent estimation and adaptive filtration]. Moscow: Nauka, 1984, 288 p.

12. L'yung L. Identifikatsiya sistem. Teoriya dlya pol’zovatelya: per. s angl. [System identification. Theory for users: translation from English]. Moscow: Nauka, 1991, 432 p.

13. Saridis Dzh. Samoorganizuyushchiesya stokhasticheskie sistemy upravleniya: per. s angl. [Self-organizing stochastic systems of control: translation from English]. Moscow: Nauka, 1980, 400 p.

14. Mosca E. Optimal, predictive and adaptive control. New Jersey: Prentice-Hall, Inc., 1995, 477 p.

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

16. Park P., Kailath T. IEEE Trans. Automat. Control. 1995, vol. 40, no. 5, pp. 895-899.

17. Jover J. M., Kailath T. Automatica. 1986, vol. 22, no. 1, pp. 43-57.

18. Hotop H.-J. Parallel Computing. 1989, no. 12, pp. 233-247.

19. Gentleman W. M. J. Inst. Math. Appl. 1973, no. 12, pp. 329-336.

20. Dyer P., McReynolds S. J. Optim. Theory Appl. 1969, no. 3, pp. 444-459.

21. Verhaegen M., P. Van Dooren IEEE Trans. Automat. Contr. 1986, vol. 31, no. 10, pp. 907-917.

22. Tsyganova Yu. V. Avtomatika i telemekhanika [Automatc and telemechanics]. 2011, no. 9, pp. 142-160.

23. Tsyganova Yu. V., Kulikova M. V. Avtomatika i telemekhanika [Automatc and telemechanics]. 2012, no. 6, pp. 34-51.

24. Tsyganova J. V., Kulikova M. V. IEEE Transactions on Automatic Control. 2013, vol. 58, no. 11, pp. 2938-2944.

25. Bjorck A. BIT. 1967, vol. 7, pp. 1-21.

26. Jordan T. L. Math. Comp. 1968, vol. 22, pp. 579-588.

Цыганова Юлия Владимировна

кандидат физико-математических наук, доцент, кафедра информационных технологий, Ульяновский государственный университет (Россия, Ульяновск, ул. Льва Толстого, 42)

E-mail: [email protected]

Tsyganova Yuliya Vladimirovna Candidate of physical and mathematical sciences, associate professor, sub-department of information technologies, Ulyanovsk State University (42 Lva Tolstogo street, Ulyanovsk, Russia)

УДК [519.254+519.654]:629.5.05 Цыганова, Ю. В.

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