авторефераты диссертаций БЕСПЛАТНАЯ БИБЛИОТЕКА РОССИИ

КОНФЕРЕНЦИИ, КНИГИ, ПОСОБИЯ, НАУЧНЫЕ ИЗДАНИЯ

<< ГЛАВНАЯ
АГРОИНЖЕНЕРИЯ
АСТРОНОМИЯ
БЕЗОПАСНОСТЬ
БИОЛОГИЯ
ЗЕМЛЯ
ИНФОРМАТИКА
ИСКУССТВОВЕДЕНИЕ
ИСТОРИЯ
КУЛЬТУРОЛОГИЯ
МАШИНОСТРОЕНИЕ
МЕДИЦИНА
МЕТАЛЛУРГИЯ
МЕХАНИКА
ПЕДАГОГИКА
ПОЛИТИКА
ПРИБОРОСТРОЕНИЕ
ПРОДОВОЛЬСТВИЕ
ПСИХОЛОГИЯ
РАДИОТЕХНИКА
СЕЛЬСКОЕ ХОЗЯЙСТВО
СОЦИОЛОГИЯ
СТРОИТЕЛЬСТВО
ТЕХНИЧЕСКИЕ НАУКИ
ТРАНСПОРТ
ФАРМАЦЕВТИКА
ФИЗИКА
ФИЗИОЛОГИЯ
ФИЛОЛОГИЯ
ФИЛОСОФИЯ
ХИМИЯ
ЭКОНОМИКА
ЭЛЕКТРОТЕХНИКА
ЭНЕРГЕТИКА
ЮРИСПРУДЕНЦИЯ
ЯЗЫКОЗНАНИЕ
РАЗНОЕ
КОНТАКТЫ


Pages:     | 1 |   ...   | 4 | 5 || 7 | 8 |

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

-- [ Страница 6 ] --

Казалось бы, основные идеи исчерпаны и все полезные результаты в этой области получены. Однако широта приложений заставила использо вать новую вычислительную среду реализации этих алгоритмов парал лельные компьютеры [8], в результате чего они оказались представлены в новой, блочной форме с ортогонализацией [17] (эти инновации изложены ниже в разд. 14). Появилось учебное пособие, рассчитанное на широкий круг пользователей и отражающее почти все достижения в этой области к концу 20 столетия [16].

Указанные исследования велись за пределами СССР. По настоящее время они остаются в России малоизвестны, применены в немногих сложных проектах (например, в авиационном приборостроении [80] или судострое нии [78, 79]) и преподаются в немногих специальных курсах [77]. В рос сийской научной литературе квадратно-корневые алгоритмы фильтрации совсем не изложены [27] или освещены слабо [64]. Так, в книге Огаркова, 1990 г. [64] этой теме посвящен лишь один параграф 3.7. Нехватка литера туры сказывается в том, что многие отечественные специалисты, работаю щие в прикладных областях регрессионного моделирования или экономет рики, продолжают использовать алгоритмы, которые можно считать уста ревшими, несмотря на то, что они испытывают значительные трудности в случае плохо обусловленной схемы наблюдения или при мультиколлине арности регрессоров [26]. Даже в публикациях последних лет [33] авторы пользуются уравнениями фильтра Калмана в их классическом виде, а не в форме современных вычислительно эффективных алгоритмов. Ситуация в России такова, что здесь лишь немногие специалисты уделяли внимание квадратно-корневым, а также ортогонализованным блочным алгоритмам, принадлежащим новому классу параллельно-ориентированных реализаций [64, 73, 78, 79, 86, 107, 108, 132, 135, 138].

13 Устойчивые алгоритмы фильтрации Первая цель данного раздела привлечь внимание студентов, маги странтов, аспирантов и специалистов к этому направлению исследований.

Здесь мы ограничиваемся лишь кратким обзором алгоритмов нового класса, поскольку более подробные обзоры доступны. Коллективная монография [73] и книга [15] содержат подробный анализ и сравнительные таблицы эффективности стандартных и квадратно-корневых алгоритмов дискретной начиная от зарождения идей у Галилео фильтрации. Прекрасный обзор, Галилея (1564–1642), через работы Карла Фридриха Гаусса (1777–1855), Нор берта Винера (1894–1964) и многих других (список огромен) до Рудольфа Эмиля Калмана (1930–) и последних достижений, датируемых до 2001 года, содержится в учебном пособии [16].

Вторая цель данного раздела показать прикладные возможности квадратно-корневых алгоритмов. Для примера рассмотрено их использова ние в актуальной задаче обнаружения момента вхождения судна в маневр, новое решение которой содержится в [41]. Соответствующий материал с опо рой на достижения в этой области [60, 90, 95, 114, 128] помещен в под разд. 13.12. Работа [95], решающая ту же задачу сопровождения маневриру ющих целей, также использует квадратно-корневую реализацию фильтров.

В отличие от [95], в [41] применен расширенный фильтр первого, а не вто рого порядка;

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

13.2 Стандартный фильтр Калмана Рассмотрим источник данных, представимый в виде линейной динами ческой системы, возмущаемой дискретным белым шумом:

x(ti+1) = (ti+1, ti)x(ti) + B(ti )u(ti) + (ti )w(ti), (13.1) z(ti ) = H(ti)x(ti) + v(ti), (13.2) где (ti+1, ti ), B(ti), (ti ), H(ti ) известные матрицы-параметры системы;

x(ti) n-мерный вектор состояния системы, u(ti) r-мерный вектор вход ного воздействия, z(ti ) m-мерный вектор измерений, w(ti ) и v(ti) неза висимые нормально распределенные векторы шумов с нулевым средним и известными ковариационными матрицами Q(ti ) и R(ti ) соответственно, при чём Q(ti) 0, R(ti) 0. Начальный вектор состояния системы x0 распре делен по нормальному закону с математическим ожиданием x0 и ковариа ционной матрицей P0.

При заданных матрицах-параметрах системы (ti+1, ti ), B(ti ), (ti ), H(ti ) 13.2 Стандартный фильтр Калмана и начальных условиях x0, P0 решение задачи оптимального оценивания вектора состояния x(ti) системы (13.1), (13.2) дается фильтром Калмана.

Алгоритм этой рекуррентной обработки данных {z(ti );

i 1} определен следующими уравнениями:

x(t) = (ti, ti1)x(t+ ) + B(ti1)u(ti1), i i1 P (ti ) = (ti, ti1)P (ti1) (ti, ti1) + (ti1)Q(ti) (ti1), + T T T T K(ti) = P (ti )H (ti )[H(ti)P (ti )H (ti ) + R(ti )], (13.3) P (t+ ) = P (t ) K(ti)H(ti)P (t ), i i i (ti) = z(ti ) H(ti )x(t), i + x(ti ) = x(ti ) + K(ti)(ti).

Замечание 13.1. Уравнения (13.1), (13.2) описывают широкий класс моделей систем, например, систем с управлением. Это могут быть и системы с возможными параметрическими нарушениями, если уравнения (13.1), (13.2) рассматривать для каждого режима функционирования (с наруше нием или нет) отдельно. Это могут быть и адаптивные системы с иденти фикацией, если матрицы, входящие в (13.1), (13.2), содержат неизвестные параметры, которые подлежат определению в процессе функционирования системы для целей слежения за состоянием системы или для целей управ ления. В частности, управление u(ti) может отсутствовать. Его наличие или отсутствие, так же как и наличие или отсутствие нарушений, не влияет на принципы эффективной численной реализации фильтров. Чтобы показать универсальность таких реализаций, управление не исключено из уравнений (13.1), (13.2). Наличие или отсутствие нарушений здесь не отмечено индек сом режима возле обозначений матриц исключительно для упрощения за писей.

Общность исходной модели (13.1), (13.2) и, соответственно, уравнений фильтра Калмана (13.3), объясняется также следующим.

Замечание 13.2. Если (ti+1, ti ) = I, u(ti) 0 и w(ti) 0, то модель (13.1), (13.2) совпадает с той моделью, которая принята в теории регрессионного моделирования, когда x(ti) = const параметр, подлежа щий оцениванию по наблюдениям z(ti ), i = 1, 2,..., N. Алгоритм (13.3) осу ществляет оптимальное оценивание этой регрессионной модели. Здесь прин ципиальная особенность заключается в том, что регрессионная модель и процесс ее оценивания последовательные. Полная матрица регрессоров H = [H(ti);

i = 1, 2,..., N ]T расщеплена на отдельные порции H(ti ) соот ветственно тем порциям наблюдений z(ti ), которые поступают в обработку 13 Устойчивые алгоритмы фильтрации в последовательные моменты времени i = 1, 2,..., N. Из теории наимень ших квадратов, изложенной выше (см. подразд. 11.6), известно, что для зна чения текущей оптимальной оценки по экспериментальным данным безраз лично, поступили эти данные в обработку целиком или порциями. Последо вательная (рекуррентная) форма алгоритма метода наименьших квадратов в форме (13.3) предпочтительнее, чем единовременная обработка всех дан ных целиком методом решения нормальных уравнений (см. разд. 12), хотя последнее все еще (к сожалению) остается традицией в эконометрике или при обработке данных наблюдений в астрономии или спутниковой гео дезии [26].

Замечание 13.3. Еще раз отметим важное обобщение, присутствую щее в модели (13.1), (13.2) и также в алгоритме фильтра (13.3). В отличие от предыдущих разделов этой книги, здесь источник обрабатываемых данных {z(ti );

i 1} не статический, а динамический. Он содержит информа цию не о неизменном по своему значению неизвестном (возможно, случай ном) векторе x(ti) = const, эта ситуация выделена в Замечании 13.2, а об изменяющемся векторе x(ti) = var состояния динамического объекта (13.1), находящегося (в общем случае) под воздействием детерминистского управления u(ti) и/или случайного возмущения w(ti).

13.3 Скаляризованная форма фильтра Калмана Когда матрица R(ti ) диагональная, R(ti) = diag [r1(ti ), r2(ti ),..., rm(ti )], возможно произвольное расщепление системы наблюдений (вектора z(ti )) на априорную и текущую части, см. выше подразд. 10.8, 11.6 или [15]. Тогда предпочтительно вместо (13.3) использовать фильтр Калмана со скалярной обработкой измерений. Скаляризуя обработку измерений в (13.3), получаем:

I. Экстраполяция:

x(t ) = (ti, ti1)x(t+ ) + B(ti1)u(ti1), i i P (t ) = (ti, ti1)P (t+ )T (ti, ti1) + (ti1)Q(ti1)T (ti1), (13.4) i i x(t+ ) = x0, P (t+ ) = P0.

0 II. Обработка измерения:

А. Начальное присваивание: P = P (t ), x = x(t ).

i i 13.4 Стабилизованный фильтр Калмана–Джозефа Б. m-кратное повторение процедуры скалярного обновления.

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

:= hT P h + rj (ti ), K := P h/, P := P KhT P, (13.5) x := x + K(z hT x) (13.6) с экстраполяцией между повторениями: P := P, x := x.

В. Завершающее присваивание: P (t+ ) := P, x(t+ ) := x, i i T где h j-й столбец матрицы H (ti );

z j-й элемент вектора z(ti ), j = 1, 2,..., m.

13.4 Стабилизованный фильтр Калмана–Джозефа Этап экстраполяции стабилизированного алгоритма совпадает с этапом экстраполяции стандартного алгоритма Калмана. Поэтому приведем лишь этап обработки измерения. Джозеф [15] предложил для алгоритма (13.3) использовать общую формулу, справедливую для матрицы P при любом, не обязательно оптимальном K: P = (I KH)P (I KH)T + KRK T. (Смысл применяемых здесь обозначений матриц виден из подразд. 13.3.) При оптимальном Kopt = P H T (H P H T + R) Замечание 13.4.

эта формула превращается в P = P Kopt H P. При подстановке сюда фор мулы для Kopt получаем дискретное алгебраическое уравнение Риккати (Discrete Algebraic Riccati Equation, DARE) [109]. Поэтому этот алгоритм может также рассматриваться как процесс решения DARE.

При таких вычислениях результирующая матрица сохраняет симметрич ность;

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

II. Обработка измерения:

А. Начальное присваивание: P = P (t ), x = x(t ).

i i Б. m-кратное повторение процедуры скалярного обновления.

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

:= hT P h + rj (ti ), v := P h, K := v/, P := P Kv T, v := P h, P := P vK T + KK T, x := x + K(z hT x) с экстраполяцией между повторениями: P := P, x := x.

13 Устойчивые алгоритмы фильтрации В. Завершающее присваивание: P (t+ ) := P, x(t+ ) := x, i i T где h j-й столбец матрицы H (ti) z j-й элемент вектора z(ti ), j = 1, 2,..., m.

13.5 Квадратно-корневой фильтр Поттера Как отмечалось в подразд. 11.8, здесь вместо матриц P (t± ), по своей при i роде положительно определенных, оперируют с их квадратными корнями S(t± ), отвечающими равенствам S(t±)S T (t±) = P (t± ). Ясно, что эти соот i i i i ношения не дают однозначного определения квадратных корней. Однако это обстоятельство в общем случае не вызывает беспокойства, поскольку эти корни однозначно определяются по разложению Холесского (например, LLT -разложение или другие, см. разд. 6).

Основная идея метода фильтрации с использованием квадратного корня состоит в замене уравнений алгоритма Калмана на аналогичные, предназна ченные для последовательного расчета S(t± ). Такой подход оправдывается i тем, что произведение S(t±)S T (t±) не теряет свойство положительной опре i i деленности (при условии полноты ранга) даже с учетом ошибок округления, тогда как ошибки округления могут приводить к потере этого свойства для матрицы P (t+ ), если она вычисляется по стандартному алгоритму (13.3).

На этапе экстраполяции, с учетом разложения матрицы P (t ) и также i Q(ti1) = Q1/2(ti1)QT /2(ti1), уравнение для нее принимает следующий вид:

S(t)S T (t) = (ti, ti1)S(t+ )S T (t+ )T (ti, ti1)+ i i i1 i + (ti1)Q1/2(ti1)QT /2(ti1)T (ti1).

Непосредственно вычислять матрицу S(t ) можно путем построения орто i гональной матрицы T размера (n + s) (n + s), такой что S T (t ) S T (t+ )T (ti, ti1) i i =T.

QT /2(ti1)T (ti1) Это обеспечивает любой алгоритм ортогональных преобразований. Хороший результат дает [80] модифицированный алгоритм Грама–Шмидта. В резуль тате после этапа экстраполяции матрица S(t) всегда будет получаться тре i угольной.

На этапе обработки измерения требуем иметь уравнение P := P KhP в виде S S T := S(In f f T /)S T. Определим число так, чтобы обеспечить тождественность: In f f T / = (In f f T )(In f f T ).

13.5 Квадратно-корневой фильтр Поттера Из получающегося отсюда квадратного уравнения, с учетом диагонально сти R(ti ) = diag [r1(ti ), r2(ti),..., rm(ti )] и равенства = f T f + rj (ti ), выби раем одно решение = (1/)/(1 + 1/), защищенное от операции вычитания положительной величины 1/ в зна менателе.

С учетом вышеизложенного получим алгоритм фильтра Поттера:

I. Экстраполяция:

x(t) = (ti, ti+1)x(t+ );

x(t+ ) = x0;

i i S T (t ) S T (t+ )T (ti, ti1) i i =T.

QT /2(ti1)T (ti1) II. Обработка измерения:

А. Начальное присваивание: S = S(t);

x = x( t ).

i i Б. m-кратное повторение процедуры скалярного обновления.

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

f := S T h;

:= f T f + rj (ti );

:= 1/(1 + 1/);

K := Sf /;

S := S Kf T ;

x := x + K(z hT x) с экстраполяцией между повторениями S := S;

x := x.

В. Завершающее присваивание: S(t+ ) := S;

x(t+) := x, i i где h j-й столбец матрицы H T (ti );

z j-й элемент вектора z(ti ), j = 1, 2,..., m.

Теперь вычисление S равносильно счету с двойной точностью P в стан дартном алгоритме Калмана. Кроме того, устранена опасность утраты мат рицей P = S S T свойства положительной определенности, что влекло бы расходимость оценок вектора состояния. Недостатком данного алгоритма является для каждого момента ti m-кратное наличие операции извлечения квадратного корня 1/.

13 Устойчивые алгоритмы фильтрации 13.6 Одноранговое обновление ПО-матриц Рассмотрим теорему об одноранговом обновлении положительно опреде ленной матрицы [15], которую доказали Agee и Turner (1972) для версии U DU T разложения ПО–матрицы. Переформулируем ее для варианта LDLT разложения.

Теорема 13.1 (Agee–Turner, 1972). Пусть матрица L нижняя тре угольная с единичной диагональю, матрица D = diag [d1, d2,..., dn ] диа T гональная и P = LDL. Пусть также заданы некоторый скаляр c и вектор– столбец a = [a1, a2,..., an ]T, такие что 0 P = P + caaT. Тогда разложение P = LDLT, где L нижняя треугольная с единичной диагональю матрица и D = diag [d1, d2,..., dn], существует и дается следующим алгоритмом:

А. Начальное присваивание: c1 = c.

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

1) di = di + ci a2 ;

i 2) Для k = i + 1, i + 2,..., n выполнять:

i. ak := ak ai lki ;

ii. ki = lki + ci ai ak /di;

A В матрицах L нетривиальные элементы l находятся только ниже диагонали.

3) ci+1 = ci di /di.

В. Завершающее присваивание: dn = dn + cn a2.

n Доказательство. Запишем квадратичную форму 0 xT P x в виде суммы полных квадратов с подстановкой в нее левой и правой частей равен ства P = P + caaT. Продолжим доказательство в виде несложного упраж нения, последовательно выделяя эти полные квадраты с правой стороны полученного равенства квадратичных форм (см. [78]). Проделайте это само стоятельно в виде упражнения. Упражнение 13.1. Аналогично предыдущему, докажите следующую версию Теоремы 13.1.

нижняя треугольная матрица и P = LLT.

Теорема 13.2. Пусть L Скаляр c и вектор–столбец a = [a1, a2,..., an ]T такие что 0 P = P + caaT.

Тогда разложение P = LLT существует и дается следующим алгоритмом:

13.7 Факторизованный фильтр Бирмана А. Начальное присваивание: c1 = c.

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

1) ii = l 2 + ci a2 ;

l ii i 2) Для k = i + 1, i + 2,..., n выполнять:

i. ak := ak ai lki/lii ;

ii. ki = lki /lii ii + ci ai /ii ak ;

l l l 3) ci+1 = ci l 2 /2.

l ii ii Упражнение 13.2.

20 c = 1.

L=, a=, 1 3 30 3 Найдите P = LLT + caaT =, L=. Алгоритм тео 06 0 ремы 13.2 дает этот же результат. Проверьте. Измените исходные дан ные, например, возьмите c = 1. Проверьте, что алгоритм дает правильный 5 результат: L =.

4/ 5 3 6/ Упражнение 13.3. Сформулируйте и докажите еще две версии теоремы об одноранговом обновлении подобно двум предыдущим теоремам, опираясь на другие разложения Холесского: P = U U T и P = U DU T [15].

13.7 Факторизованный фильтр Бирмана Основная идея этого алгоритма состоит в разложении ковариацион ной матрицы P в произведение двух треугольных матриц и диагональ ной матрицы между ними. Можно рассматривать два варианта алго ритма: LD-алгоритм или U D-алгоритм Бирмана. Последний изложен в [15]. Здесь рассмотрим LD-алгоритм, т. е. используем разложение P (t± ) = i = L(t± )D(t±)LT (t± ), где L(t±) нижние треугольные матрицы с единичной i i i i ± диагональю, D(ti ) диагональные матрицы.

13 Устойчивые алгоритмы фильтрации На этапе экстраполяции с учетом разложения матрицы Q(ti1) = = Lq (ti1)Dq (ti1)LT (ti1) и аналогичного разложения матрицы P (t ) урав q i нение для P (ti ) принимает следующий вид:

P (t ) = L(t)D(t)LT (t) = (ti, ti1)L(t+ )D(t+ )LT (t+ )T (ti, ti1) + i i i i i1 i1 i + (ti1)Lq (ti1)Dq (ti1)LT (ti1)T (ti1).

q Представим матрицу P (t ) в следующем виде:

i P (t ) = (ti, ti1)L(t+ ) (ti1)Lq (ti1) i i LT (t+ )T (ti, ti1) Diag [D(t+ ), Dq (ti1)] i, Lq (ti1)T (ti1) T i т. е. P (t ) = W (ti )D(ti)W T (ti), i T w1 (ti ) =..., (ti, ti1)L(t+ ) W (ti) = (ti1)Lq (ti1) i T wn (ti ) D(ti ) = Diag [D(t+ ), Dq (ti1)] = Diag [D1 (ti),..., DN (ti )], i где W (ti) – матрица размера (n(n+s)) и N = n+s, s размерность шума w(ti ) в уравнении состояния (13.1). Факторы L(t ) и D(t ) вычисляются i1 i методом модифицированного взвешенного алгоритма Грама–Шмидта [15] на этапе экстраполяции.

Перефразируя его для нижних треугольных факторов, запишем следу ющий результат (для удобства изложения будем писать L, L, D, D, D, W вместо L(t), L(t+), D(t ), D(t+ ), D(ti ), W (ti) соответственно):

i i i i Теорема 13.3 (Взвешенная ортогонализация Грама–Шмидта и фак торизация матрицы аналог теоремы VI.4.1 из [15]). Пусть векторы N w1,..., wn образуют линейно независимую систему, wi R, N n и D = = Diag [D1,..., Dn] 0. Применим следующий алгоритм (нижний индекс указывает номер нетривиального элемента матрицы / вектора):

(1) А. Начальное присваивание: vk = wk, k = 1,..., n.

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

(j1) (j1) dj1 = (vj1 )T Dvj1 ;

(j1) T (j1) lk,j1 = (vk ) Dvj1 /dj1, k = j,..., n;

(j) (j1) (j1) lk,j1vj1, vk = vk k = j,..., n.

13.7 Факторизованный фильтр Бирмана (n) (n) В. Завершающее присваивание: dn = (vn )T Dvn.

(j) Тогда vj = vj, j = 1,..., n, где vj взвешенные взаимно ортогональные T векторы: vi Dvj = dj i,j (i,j символ Кронекера) и T T w1 v... D w1... wn = L... D v 1... v n LT = LD LT.

T T wn vn Из этого алгоритма имеем результат: L(t) = L, D(t ) = D.

i i Этап обработки измерений представим в виде теоремы, формулируя ее для нижних треугольных сомножителей (LD-версия) [78] вместо U D-версии теоремы Бирмана из [15].

Теорема 13.4 (Алгоритм Бирмана). Пусть калмановская процедура скалярного обновления (13.5), (13.6) использует разложения P = LDLT и P = LDLT, где L нижние треугольные с единичной диагональю матрицы, D диагональные (с положительными элементами) матрицы. Тогда данная процедура (13.5), (13.6) эквивалентна пунктам В, Г, Д и Е, следующего алгоритма.

II. Обработка измерения:

А. Начальное присваивание: L = L(t ), D = D(t ), x = x(t).

i i i Б. m-кратное повторение процедуры скалярного обновления.

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

В. Вычислить векторы f = [f1, f2,..., fn]T = LT h;

v = [v1, v2,..., vn]T = Df.

Г. Задать начальные значения = r;

K = [0... 0 vn ]T.

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

начало := + vifi ;

:= 1/ ;

:= fi ;

di := di ;

li := li + K ;

K := K + li vi ;

(A) :=.

конец Е. Вычислить векторы := (z hT x) ;

x := x + K с экстраполяцией между повторениями:

L := L ;

D := D ;

x := x.

13 Устойчивые алгоритмы фильтрации Ж. Завершающее присваивание:

L(t+) := L ;

D(t+ ) := D ;

x(t+ ) := x.

i i i j-й столбец матрицы H T (ti );

z Здесь h j-й элемент вектора z(ti );

r j-й элемент rj (ti ) диагональной матрицы ковариаций шума измерений R(ti), j = 1, 2,..., m номер скалярного изме рения в составе вектора измерений z(ti ) в момент ti.

Замечание 13.5. Строка A в алгоритме на стр. 281 выглядит так:

Для k = i + 1 до n выполнять начало lki := lki + Kk, Kk := Kk + lkivi конец Она пропускается при i = n. Здесь Kk есть k-й элемент того вектора K, который существует в цикле Д алгоритма на стр. 281.

Доказательство. Применим для (13.5) факторизацию Холесского вида P = LDP T и P = LD LT. Получим LDP T = L(D + cvv T )LT с обозначениями f = LT h, c = 1/, v = Df.

Применяя теорему 13.1 к выражению D + cvv T, в ней следует считать, что:

P D, c 1/ и a v. Отсюда получаем LD LT = LDLT + cvv T при том, что в этом выражении имеем L = I. Когда L и D будут найдены, мы получим LDP T = LLD LT LT = D = D, L = LL. (13.7) Запишем теорему 13.1 применительно к выражению D + cvv T. Получим:

А. Начальное присваивание: c1 = c.

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

1) di = di + ci vi. x 2) Для k = i + 1, i + 2,..., n выполнять:

i. vk := vk vilki = vk. r Здесь L = I, т. е. lki = 0. Это действие исчезает.

ii. ki = lki + (civi /di)vk = (ci vi/di)vk. y l 13.7 Факторизованный фильтр Бирмана 3) ci+1 = ci di /di. z { В. Завершающее присваивание: dn = dn + cn vn.

Этот алгоритм теоретически верен, но численно неустойчив из-за того, что должно выполняться условие: i : ci 0. В силу этого величины di в x могут стать отрицательны, что противоречит требованию положитель ной определенности P 0. Устраним этот недостаток, эквивалентно пре образуя алгоритм на стр. 282–283 по пунктам x, y, z, {. Преобразуем x:

2 di /ci = di /ci + vi. С учетом z: di/ci+1 = di/ci, поэтому di/ci+1 = di/ci + vi, следовательно, 1/ci+1 = 1/ci + vi /di. Однако vi = di fi. Поэтому 1/ci = = 1/ci+1 + vi fi. Введем обозначение: i : i = 1/ci. Это позволяет запи сать:

i = i+1 + vifi, i = n 1, n 2,..., 2, 1. (13.8) Финальное значение, т. е. 1, согласно введенному обозначению, должно сов пасть с = r + v T f, где r = r(tj ), так как c1 = c и c = 1/. Чтобы это произошло, в алгоритме необходимо стартовать от значения n = r + vnfn.

Таким образом, из x получили алгоритм (13.8), заменяющий z.

Теперь выведем из z алгоритм, заменяющий x:

i = n 1, n 2,..., 2, 1.

di = di (ci /ci+1) = di (i+1/i ), Осталось преобразовать y. Из y имеем ki = i vk, i = (ci vi/di), а из l z ci /di = ci+1/di, то есть i = fi/i+1, так как vi = di fi. Следовательно, алгоритм для ki вместо y приобретает вид:

l Для i = n 1, n 2,..., 2, 1 вычислять i = fi /i+1.

Для k = i + 1 до n вычислять ki = i vk.

l Наконец, преобразуем последнее действие { на стр. 283. Последовательно находим:

2 dn /cn = dn /cn + vn, (dn /dn)(1/cn) = (1/cn) + (vn/dn ), (vn/dn ) = (dn fn)2 /dn = vn fn.

Из обозначения (1/ci+1) = i+1 при i = n 1 следует (1/cn) = n, т. е.

n = (dn /dn)(1/cn) + vn fn. С другой стороны, известно стартовое значе ние n = r + vnfn, следовательно, dn = dn r/n. Запишем результирующий эквивалентный алгоритм:

13 Устойчивые алгоритмы фильтрации Задать начальные значения: n = r + vn fn, dn = dnr/n.

Для i = n 1, n 2,..., 2, 1 вычислять i = i+1 + vi fi ;

di = di (i+1/i ) ;

i = fi/i+1.

Для k = i + 1 до n вычислять ki = i vk.

l Отмеченный выше недостаток устранен: теперь i : di 0. Проведен ное эквивалентное преобразование привело к инверсии направления вычис лений. В исходном массиве для нетривиальных элементов матриц L и D исходные столбцы сначала имеют (начиная с диагонали вниз) вид:

d1 d2 dn,,...,,.

dn l1 l2 ln После выполнения этого алгоритма они имеют следующее наполнение:

d1 d2 dn,,...,,.

dn 1 2 n l l l Вернемся к последней формуле в (13.7), чтобы в согласии с ней перейти от промежуточной матрицы L к окончательной матрице L.

Матрица L нижняя треугольная с единицами на главной диагонали.

В ней поддиагональная часть i-го столбца, согласно нижней строке алго ритма на стр. 284, найдена в виде iT = i [vi+1, vi+2,..., vn], где использо l T ваны элементы из вектора v = [v1, v2,..., vn], начиная от элемента vi+1 и далее. Введем вспомогательные обозначения n-мерных векторов:

0 0 0 v.. v v..

.. 2 (n1) v (1) v (n) (2) v = 0, v = 0,..., v = 3, v = 3 = v.

..

..

0 vn1..

vn vn vn vn Теперь видно, что L = I + 1 v (2) 2 v (3) · · · n1 v (n) 0.

13.8 Квадратно-корневой фильтр Карлсона Отсюда L = LL = L + 1 Lv (2) 2 Lv (3) · · · n1Lv (n) 0.

Введем обозначения возникающих здесь матриц-столбцов:

K2 = Lv (2), K3 = Lv (3), · · ·, Kn = Lv (n).

В качестве начального столбца имеем Kn = Lv (n). Далее Kn1 = Lv (n1) = Kn + ln1vn1, Kn2 = Lv (n2) = Kn1 + ln2vn2, ···, Ki = Ki+1 + li vi, i = n 1, n 2,..., 2, 1, и в конце получаем K1 = Lv (1) = Lv = K, K вспомогательное обозначение, причем здесь для строгости записей надо понимать, что li, i = n1, n2,..., 2, 1, обозначает весь j-й столбец матрицы L вместе с тривиальными элементами 0 и 1. Данный в теореме 13.4 алгоритм не содержит операции извлечения квад ратного корня, а работа с треугольными матрицами требует меньшего числа арифметических операций по сравнению с обычными.

13.8 Квадратно-корневой фильтр Карлсона Этап экстраполяции здесь в точности совпадает с этим этапом в алго ритме Поттера. Выведем только алгоритм этапа обработки измерения.

Упражнение 13.4. Выведите алгоритм Карлсона обработки измере ния. Вывод можно сделать в полном соответствии с тем, как выше была доказана теорема 13.4 (см. стр. 281). Другой способ вывода как следствие теоремы 13.4: за L принять формально то, что там было LD1/2. Первый спо соб предпочтительнее для понимания доказательства, второй полезен для освоения техники вычислений. Получите следующий результат.

Теорема 13.5. Пусть калмановская процедура скалярного обновле ния (13.5), (13.6) использует разложения P = LLT и P = LLT, где L и L нижние треугольные матрицы. Тогда данная процедура (13.5), (13.6) эквивалентна пунктам В, Г, Д и Е, следующего алгоритма.

13 Устойчивые алгоритмы фильтрации II. Обработка измерения:

А. Начальное присваивание: L = L(t ), x = x(t).

i i Б. m-кратное повторение процедуры скалярного обновления.

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

В. Вычислить векторы f = [f1, f2,..., fn]T = LT h.

Г. Задать начальные значения: = r;

K = [0... 0 lnn fn ]T.

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

начало := + fi2 ;

:= / ;

:= fi/() ;

lii := lii ;

li := li + K ;

K := K + li fi ;

(A) :=.

конец Е. Вычислить векторы := (z hT x)/ ;

x := x + K с экстраполяцией между повторениями:

L := L ;

x := x.

Ж. Завершающее присваивание:

L(t+) := L ;

x(t+ ) := x.

i i j-й столбец матрицы H T (ti );

z Здесь h j-й элемент вектора z(ti );

r j-й элемент rj (ti ) диагональной матрицы ковариаций шума измерений R(ti), j = 1, 2,..., m номер скалярного изме рения в составе вектора измерений z(ti ) в момент ti.

Замечание 13.6. Строка A в алгоритме на стр. 286 выглядит так:

Для k = i + 1 до n выполнять начало lki := lki + Kk, Kk := Kk + lkifi конец Она пропускается при i = n. Здесь Kk есть k-й элемент того вектора K, который существует в цикле Д алгоритма Карлсона на стр. 286.

13.9 Редуцированный фильтр Бирмана 13.9 Редуцированный фильтр Бирмана В приложениях выбор измерительных средств бывает ограничен так, что в измерение z попадает лишь часть (заранее известная) элементов оценива емого вектора x. Пусть эта часть первые qj n элементов для j-й строки матрицы наблюдений:

hT = [ ··· 0···0] = [ ··· 0 · · · 0 ]. (13.9) j qj n, hqj = Теорема 13.6 ([78]). Пусть калмановская процедура скалярного обновления (13.5), (13.6) использует разложения P = LD LT и P = LDLT, где L нижние треугольные с единичной диагональю матрицы, D диаго нальные (с положительными элементами) матрицы. Тогда данная процедура (13.5), (13.6) эквивалентна пунктам В, Г, Д и Е, следующего алгоритма.

II. Обработка измерения:

А. Начальное присваивание: L = L(t ), D = D(t ), x = x(t).

i i i Б. m-кратное повторение процедуры скалярного обновления.

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

В. Вычислить векторы f = LT h ;

v = Df. Их вид:

0 · · · 0 ]T, v = [ 0 · · · 0 ]T.

··· ··· f =[ Г. Задать начальные значения = r;

K = [0... 0]T.

Д. Для i = q + 1, q + 2,..., n выполнять: Ki := liq vq.

Е. Для i = q, q 1,..., 2, 1 выполнять:

начало := + vifi ;

di := di / ;

:= fi/ ;

Ki := vi ;

li := li + K ;

K := K + livi ;

(A) :=.

конец Ж. Вычислить векторы := (z hT x)/ ;

x := x + K с экстраполяцией между повторениями:

L := L ;

D := D ;

x := x.

З. Завершающее присваивание:

L(t+) := L ;

D(t+ ) := D ;

x(t+ ) := x.

i i i 13 Устойчивые алгоритмы фильтрации Здесь hT j-я строка матрицы H(ti ), имеющая вид (13.9);

z j-й элемент вектора z(ti );

r j-й элемент rj (ti ) диагональной мат рицы ковариаций шума измерений R(ti), j = 1, 2,..., m номер скалярного измерения в составе вектора измерений z(ti ) в момент времени ti.

Замечание 13.7. Строка A в алгоритме на стр. 287 выглядит так же, как в замечании 13.5 на стр. 282.

Замечание 13.8. Представление (13.9) мотивирует использование LD-версии фильтра Бирмана по подразд. 13.7. Если же hT = [ 0 · · · 0 · · · ], j то сокращенный объем вычислений будет получаться при использовании U D-версии фильтра Бирмана [15].

Замечание 13.9. Сокращение объема вычислений в редуцирован ных таким образом версиях получается наибольшим среди алгоритмов этого класса, поскольку значение q берется как qi индивидуально для каждого из m повторений цикла Е в алгоритме на стр. 287.

Замечание 13.10. Массив D может быть общим для D и D, и в нем обновляются только первые q элементов. В массивах для матриц L и L обновляются только первые q 1 столбцов, а столбцы с q-го по n-й не обновляются и могут быть общими.

13.10 Редуцированный фильтр Бар-Ицхака Если в исходном фильтра Калмана (13.3) свойство (13.9) рассматривать сразу для всех m строк матрицы наблюдений H(ti) = H = [ H mq 0 ], (13.10) это может служить поводом для сокращения объема вычислений, как это произошло выше в подразд 13.9. В (13.10) H mq ненулевая подматрица, mq указывает ее размер (m q), где нулевая подматрица, ее размер (m s), s = n q.

Свойство (13.10) означает, что оцениваемый вектор распадается на две части: x = [ xq xs ]. Часть xq попадает в вектор измерений z = z(ti ), а часть 13.11 Редуцированный фильтр Бар-Ицхака–Медана xs нет. Соответственно этому, каждую P -матрицу, т. е. P и P, рассмотрим поблочно как матрицу следующего вида [89]:

P qq (P sq )T P=. (13.11) P sq P ss Теорема 13.7 (Bar-Itzhack, 1980). Если выполнено условие (13.10), алгоритм (13.3) распадается на независимый редуцированный фильтр раз мерности q для измеряемых компонент xq вектора x (аргумент дискретного времени ti для простоты опущен):

K = P (H ) H P (H ) + R, mq T mq T qm qq mq qq (13.12) P qq = P qq K qmH mq P qq, xq = xq + K qm (z H mq xq ) и фильтр размерности s = nq, зависимый от предыдущего фильтра (13.12), для неизмеряемых компонент xs вектора x :

sq sq qq K =P P, sq sq qq P =K P, (13.13) sq T (K ), ss ss sq qq qq P =P K P P s s sq q q x = x + K (x x ).

Доказательство. Подстановка выражений (13.10) и (13.11) в уравнения стандартного фильтра Калмана (13.3) после несложных алгебраических пре образований приводит к (13.12) и (13.13). Полученный алгоритм лишь выделяет редуцированный фильтр (13.12), но для него задача LD-факторизации остается актуальной. Ее решает сле дующий алгоритм (подразд. 13.11).

13.11 Редуцированный фильтр Бар-Ицхака–Медана Этот результат представлен в 1983 году как новый LD-алгоритм обнов ления оценок по измерениям. Он является лучшим по сравнению с обычным LD-алгоритмом в случае, когда число элементов вектора состояния, непо средственно попадающих в вектор измерений, меньше, чем размерность век тора состояния. Такая ситуация типична для аэрокосмических приложений.

13 Устойчивые алгоритмы фильтрации Так, в инерциальных навигационных системах число состояний может быть более 40, а число измеренных состояний только 3.

Теорема 13.8 (Bar-Itzhack–Medan, 1983 [89]). Пусть при условии (13.10) алгоритм Калмана (13.3) использует разложения P = LD LT и P = = LD LT с обозначениями P = P (t ) и P = P (t+ ), причем P и P рассмот i i рены поблочно, как в (13.11), и L, D, L, D разложены по типу выражений Lqq Dq 0 L=, D=. (13.14) Lsq Lss Dss Тогда этот алгоритм эквивалентен следующему алгоритму (излагаемое отно сится лишь к этапу II обработки измерения с матрицей вида (13.10)):

II. Обработка измерения, эквивалентная работе алгоритма (13.12):

А. Начальное присваивание: L = Lqq (t ), D = Dq (t ), x = xq (t).

i i i Б. m-кратное повторение процедуры скалярного обновления.

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

В. Вычислить векторы f = [f1, f2,..., fq ]T = LT h;

v = [v1, v2,..., vq ]T = Df.

Г. Задать начальные значения = r;

K = [0... 0 vq ]T.

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

начало := + vifi ;

:= 1/ ;

:= fi ;

di := di ;

li := li + K ;

K := K + li vi ;

(A) :=.

конец Е. Вычислить векторы := (z hT x) ;

x := x + K с экстраполяцией между повторениями:

L := L ;

D := D ;

x := x.

Ж. Завершающее присваивание по п. II:

Lqq (t+ ) := L ;

Dq (t+) := D ;

xq (t+ ) := x.

i i i j-й столбец подматрицы (H mq )T (ti) из (13.10);

z Здесь h j-й элемент вектора z(ti );

r j-й элемент rj (ti) диагональной матрицы 13.11 Редуцированный фильтр Бар-Ицхака–Медана ковариаций шума измерений R(ti ), j = 1, 2,..., m номер скаляр ного измерения в составе вектора измерений z(ti ) в момент ti.

III. Обновление оценок для неизмеряемых компонент xs вектора x, эквивалентное работе алгоритма (13.13):

З. Вычислить: K sq = Lsq Lqq, Lsq = K sq Lqq, (13.15) ss ss s s L =L, D =D, xs = xs + K sq (xq xq ).

И. Завершающее присваивание по п. III:

Lsq (t+) := Lsq ;

Lss (t+) := Lss ;

Ds (t+ ) := Ds ;

xs (t+) := xs.

i i i i Замечание 13.11. Строка A в алгоритме на стр. 290 выглядит так:

Для k = i + 1 до q выполнять начало lki := lki + Kk, Kk := Kk + lkivi конец Она пропускается при i = q. Здесь Kk есть k-й элемент того вектора K, который существует в цикле Д алгоритма на стр. 290.

Подстановка (13.14) в разложение P = LDLT дает Доказательство.

P qq = Lqq Dq (Lqq )T, P sq sq q qq T (13.16) = L D (L ), P ss = Lsq Dq (Lsq )T + LssDs (Lss)T.

Из первого уравнения (13.13) и из первых двух уравнений (13.16), учиты вая добавление верхней тильды · над любым символом ·, получаем первое уравнение из (13.15). Из второго уравнения (13.13) и из первых двух урав нений (13.16), учитывая добавление верхней крышки · над любым символом ·, получаем второе уравнение из (13.15), имея в виду, что P 0 влечет P qq 0, Dq 0 и det Lqq = 0. Из третьего уравнения (13.13), из первых двух уравнений (13.15) и из первого и третьего уравнений (13.16) следует Lss Ds (Lss )T = LssDs (Lss)T.

13 Устойчивые алгоритмы фильтрации Благодаря единственности LD-факторизации, это дает третью строку из выражений (13.15). Данный алгоритм привлекателен своей простотой, но применять его надо с большой осторожностью, принимая во внимание следующие замечания.

Замечание 13.12. Из первых двух уравнений (13.13) видно, что P (P qq )1 = P sq (P qq )1, так же как первые два уравнения из (13.15) дают sq Lsq (Lqq )1 = Lsq (Lqq )1. Если на этапе экстраполяции должно быть сохра нение ковариаций, что соответствует операторам P sq := P sq и P qq := P qq, то значение K sq не будет вообще изменяться. Известно, что такое сохранение ковариаций возникает в задаче чистой регрессии в предположении, что оцениваемый вектор x постоянен во времени. Поэтому данный алгоритм в задаче регрессионного моделирования неприменим. Это понятно также из того, что в этой задаче условие (13.10) означает полную ненаблюдаемость компонент xs оцениваемого вектора x = const.

Замечание 13.13. Для алгоритма Бар-Ицхака (см. подразд 13.10) из P sq = 0, а здесь из Lsq = 0, следует K sq = 0, xs = xs, а также P sq = 0, P ss = P ss и Lsq = 0. Поэтому, если в измерение z непосред ственно попадает только часть xq оцениваемого вектора x и если априори выбрана такая матрица P (для данного алгоритма L), что в ней P sq = 0 (со ответственно, Lsq = 0), то в принципе нельзя улучшить априорную оценку xs. Отсюда видно, насколько важно задавать P sq = 0 (соответственно, Lsq = 0) до запуска алгоритма. В динамических задачах оценивания это не столь критично, так как там на этапе экстраполяции вычисляют новое P := P T + QT, а не просто полагают P := P.

Замечание 13.14. Выбор LD-факторизации здесь не произво лен, как отмечается в [89], а подчинен желанию обеспечить независимость редуцированного фильтра (13.12) (соответственно, этапа II алгоритма на стр. 290), а также правым, а не левым размещением нулевого блока в (13.10).

Замечание 13.15. Первое выражение в (13.15) реализуется проще без вычисления обратной матрицы, если его переписать как линейную сис тему (Lqq )T (K sq )T = (Lsq )T и решать ее обратной подстановкой относительно (K sq )T.

Следствие 13.1. В условиях теоремы 13.8 ее утверждение справед 13.12 Задача сопровождения судна на траектории ливо для алгоритма, в котором (13.15) заменены следующими выражениями:

Lqq = (Lqq )1Lqq, sq qq sq L =L L, ss ss s s (13.17) L =L, D =D, q = (Lqq )1(xq xq ), s s sq q x =x +L.

Замечание 13.16. Первое выражение из (13.17) реализуется не вы числением обратной матрицы, а применением процедуры обратной подста новки к системе Lqq Lqq = Lqq относительно Lqq. Аналогично, четвертое выра жение из (13.17) реализуется применением обратной подстановки к системе Lqq q = xq xq относительно q. Следствие 13.1 дает более эффективный алгоритм, чем теорема 13.8, когда q + 1 s = n q, т. е. при q (n 1)/2.

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

Модель движения судна Текущее состояние судна на плоской траектории с точки зрения кинема тики (а не динамики) движения, характеризуется пятью параметрами: две координаты, скорость, направление и угловая скорость движения. Модель пространственного движения включает 12 уравнений для координат центра масс, углов Эйлера, составляющих линейной и угловой скоростей в связан ной системе координат [60]. Представим движение судна как дискретный стохастический процесс sk, подчиняющийся следующему уравнению [114]:

x + (2/) v sin(t/2) cos( + t/2) y + (2/) v sin(t/2) sin( + t/2) sk = v + wk1. (13.18) + t k 13 Устойчивые алгоритмы фильтрации Модель (13.18) есть модель состояния sk, что позволяет применять тео рию фильтрации Калмана. Вектор sk = (x, y, v,, )T содержит пять вели чин: x, y географические координаты цели (широта и долгота в прямо угольной системе координат в метрах), v линейная скорость (м/с), направление движения (курсовой угол в радианах от северного направления по часовой стрелке), угловая скорость (рад/с). Аддитивно действую щий белый шум с единичной ковариационной матрицей Q(ti) = I обозначен через wk1, а t период дискретизации времени. Ковариация шума w равна [114] t3 t T T T, t2v, 2 3 cov(w) E ww = = diag, t2 00 t 2 где v и дисперсии случайных гауссовских процессов, аддитивно воз действующих на скорость корабля и его угловую скорость, соответственно.

Модель движения (13.18) включает угловую скорость и потому пригодна для анализа траектории как маневрирующего судна, так и судна, движу щегося равномерно и прямолинейно (в последнем случае угловая скорость судна-цели равна нулю).

Модель измерений Будем считать, что наблюдения за судном ведутся при помощи радиоло кационной станции (РЛС), возвращающей измерения сферических (на плос кости полярных) координат цели: азимут z и дальность до цели z.

Построим модель измерения в декартовых координатах. Будем называть такие преобразованные измерения псевдоизмерениями. Первичные измере z ния, приходящие от РЛС, обозначим zsph =. Тогда z = +, z z = +, где и истинные значения дальности и пеленга, а и соответствующие погрешности измерения. Обозначим шум измерений полярных координат vsph =. Его ковариационная матрица имеет вид T Rsph E vsphvsph = 2, т. е. матрица диагональна, а шумы по даль ности и азимуту взаимно не коррелированы (ошибки измерения дальности не связаны с ошибками измерения направления на цель).

13.12 Задача сопровождения судна на траектории Перейдем от измерений в сферической (точнее, полярной) системе коор динат к псевдоизмерениям в декартовых координатах:

zx z cos(z ) zdec = =.

zy z sin(z ) Для удобства используем здесь и далее индекс dec в случаях работы с декартовыми координатами. Малые приращения координат в сферической и декартовой системе координат связаны соотношением:

cos(z ) z sin(z ) dzx dz dzdec = =.

dzy sin(z ) z cos(z ) dz cos(z ) z sin(z ) Обозначим T =. Перейдем от дифференциалов к sin(z ) z cos(z ) конечным разностям и получим закон преобразования случайных погреш ностей измерений из полярных координат в декартовы:

cos(z ) z sin(z ) zx z zdec = = = T zsph.

zy sin(z ) z cos(z ) z Ковариационная матрица погрешностей измерений в декартовых координа тах равна E zdec zdec = E T zsph zsph T T = T Rsph T T.

T T Rdec Матрица T здесь зависит от времени, однако для простоты изложения ин декс времени опускаем.

Выразим псевдоизмерения через истинные декартовы координаты цели и погрешности псевдоизмерений:

z cos(z ) ( + ) cos( + ) zdec = = z sin(z ) ( + ) sin( + ) cos() + zx x zx = +.

sin() + zy y zy x Здесь обозначим: s = точные значения декартовых координат цели, y zx v= – погрешности псевдоизмерений декартовых координат. Тогда zy zdec = s+v. Вследствие нелинейности преобразования координат ковариаци онная матрица погрешностей псевдоизмерений Rdec не является диагональ ной. Приведем ковариационную матрицу к единичному виду, т. е. выполним декорреляцию погрешностей псевдоизмерений.

13 Устойчивые алгоритмы фильтрации T Представим матрицу Rdec как произведение Rdec = Sdec Sdec, где S = 1/2 1 1 1 1 = T Rsph. Обозначим zdec = Sdec zdec = Sdec (s + v) = Sdec s + Sdec v = Sdec s + v.

1 1 T E v v T = E Sdec v v T (Sdec ) Найдем ковариацию v : cov() = Rdec v = 1 T = Sdec Rdec Sdec = I. Будем называть zdec нормализованными псевдоизмере ниями, а v вектором погрешностей нормализованных псевдоизмерений (его компоненты взаимно не коррелированы, так как Rdec = I ). Отсутствие корреляции компонентов вектора погрешностей позволяет применять ска ляризованную обработку измерений.

Поскольку уравнения в системе (13.18) обладают значительной нелиней ностью, непосредственное применение линейного фильтра Калмана, CKF, не представляется возможным. Одним из способов обойти это ограничение является использование расширенного фильтра Калмана (Extended Kalman Filter, EKF).

Стандартный алгоритм расширенного фильтра Калмана Расширенный фильтр Калмана является субоптимальным алгоритмом [90]. В области фильтрации и идентификации систем этот термин не обла дает уникальностью смысла. В терминологии данного подраздела расши ренность в сравнении с линейным фильтром состоит в возможности при нятия нелинейной модели движения цели и/или модели измерений. В дан ном контексте это общепринятая терминология. Формулы расширенной фильтрации получены в работе [128] методом инвариантного погружения.

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

x(ti+1) = f [x(ti)] + w(ti), z(ti ) = h[x(ti)] + v(ti).

Функция f [·] вычисляет состояние системы в момент времени ti+1 по состоянию в момент времени ti, функция h[·] преобразует вектор состояния к виду, в котором измерения поступают на вход фильтра, w(ti) и v(ti) шум процесса и шум измерения в момент ti, соответственно. Когда обе функции f [·] и h[·] линейные, расширенный фильтр превращается в стандартный фильтр Калмана.

13.12 Задача сопровождения судна на траектории Расширенный фильтр Калмана, как и стандартный, содержит два этапа:

этап 1 экстраполяция и этап 2 обработка измерения.

Экстраполяция (здесь учтено, что Q(ti) = I):

x(ti) = f [(t+ )], x i f F (ti1) =, x x=(t+ ) x i P (t ) = F (ti1)P (t+ )F T (ti1) + (ti1)T (ti1).

i i Обработка измерения (здесь учтено, что R(ti ) = Rdec = Im):

(ti) = z(ti ) h[k (t )], xi h H(ti) =, x x=(t+ ) x i K(ti) = P (t )H T (ti )[H(ti)P (t)H T (ti ) + Im]1, i i P (t+ ) = P (t ) K(ti)H(ti)P (t ), i i i x(t+) = x(t) + K(ti)(ti).

i i Приведенные уравнения расширенной фильтрации по структуре совпа дают с уравнениями обычного, линейного фильтра. Это следствие того, что в них сохранены только первые (линейные) члены разложения функ ций f [·] и h[·] в ряды Тейлора. В силу этого, рассматриваемый расширен линейный (так называемый расширенный фильтр первого ный фильтр порядка).

Линеаризация моделей процесса и измерений производится относительно последней (текущей) оценки вектора состояния. Для этого вычисляют две матрицы Якоби: F (ti1) и H(ti), которые затем подставляют в уравнения стандартной калмановской фильтрации.

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

13 Устойчивые алгоритмы фильтрации Пусть на вход фильтра поступают нормализованные псевдоизмерения (обладающие, как показано выше, единичной ковариационной матрицей ошибок). Для модели процесса, описываемой уравнениями (13.18), функции f [·] и h[·] имеют следующий вид:

x + (2/) v sin(t/2) cos( + t/2) y + (2/) v sin(t/2) sin( + t/2) f [s(ti)] =, v + t ti x h[s(ti )] =.

y ti Обозначим: = + t. Вычислим матрицы Якоби:

f F (ti1) F= x x=(t+ ) x i 1 0 sin()sin() v(cos()cos()) v(t cos()+sin()sin()) cos()cos() v(sin()sin()) v(t sin()+cos()cos()) 0 1 =0 0, 1 0 0 0 0 1 t 00 0 0 1 t+ i h H(ti ) =H.

x x=(t+ ) x i Поскольку нормализованные псевдоизмерения на входе фильтра про изводятся (т. е. формально представлены) в той же системе координат, в которой записаны координаты объекта в векторе состояния, функция h[·] выполняет линейное преобразование, и матрица H(ti ) = H не зависит от времени. Приближенность записи H(ti ) = H здесь состоит в том, что в усло виях малости величин zsph косинус угловой погрешности заменен на 1.

При прямолинейном равномерном движении ( = 0) функции f [·] и F принимают следующую форму:

13.12 Задача сопровождения судна на траектории x + (2/) v sin(t/2) cos( + t/2) y + (2/) v sin(t/2) sin( + t/2) f=0[s(ti)] =lim v = 0 + t ti x + vt cos() y + vt sin() =v, 0 ti v t2 sin() 1 0 t cos() v t sin() v t2 cos() 0 1 t sin() v t cos() F=0(ti1) = 0 0.

1 0 0 0 0 1 t 00 0 0 1 t+ i Конкретные значения элементов F во время фильтрации получаются подстановкой сюда элементов соответствующего вектора оценки состояния, однако для упрощения записей в последующих формулах символы экстра поляционной оценки (тильда ) и отфильтрованной оценки ( крышка ) опущены над всеми входящими величинами.

Учывая изложенное, скаляризованный алгоритм Поттера для нелинейной фильтрации данных о траектории судна–цели получим в следующем виде:

I. Экстраполяция:

x + (2/) v sin(t/2) cos( + t/2) y + (2/) v sin(t/2) sin( + t/2) s(t) = v, s(t+ ) = s0, i + t t+ i sin()sin() v(cos()cos()) v(t cos()+sin()sin()) 10 cos()cos() v(sin()sin()) v(t sin()+cos()cos()) 0 1 F (ti, ti1) = 0 0, 1 0 0 0 0 1 t 00 0 0 1 t+ i 13 Устойчивые алгоритмы фильтрации S T (t ) S T (t+ )F T (ti, ti1) i i =T.

T (ti ) II. Обработка измерения:

А. Начальное присваивание: S = S(t);

s = s( t ).

i i Б. m-кратное повторение процедуры скалярного обновления:

q := S T h, := q T q + 1, := 1/(1 + 1/), K := Sq/, S := S Kq T, s := s + K(z hT s) с экстраполяцией между повторениями S := S, s := s.


В. Завершающее присваивание: S(t+ ) := S, s(t+ ) := s, i i T где h j-й столбец матрицы H (ti);

z j-й элемент вектора z(ti ) (j = 1, 2,..., m).

Замечание 13.17. В рассматриваемой прикладной задаче может быть использован также расширенный фильтр второго порядка, как это сде лано в работе [95]. Важно отметить, что для придания своему фильтру чис ленной устойчивости авторы [95] также приводят его к квадратно-корневой форме.

13.13 Пример задачи с мультиколлинеарностью Пусть требуется решить несовместную систему Ax z при 1 2 11 1 21 1 2 A =..., z =.. (13.19).....

.....

N11 NN по критерию наименьших квадратов: zAx 2 = (zAx)T (zAx) minx.

Этот пример здесь выбран потому, что он иллюстрирует тяжелый слу чай мультиколлинеарности: столбцы матрицы A в (13.19) не просто колли неарны, они совпадают! Данный пример в его статистической интерпрета ции эквивалентен решению задачи оптимального оценивания неизвестного параметра x по экспериментальным данным N z = A + v, v N (0, I), xN (13.20) 13.13 Пример задачи с мультиколлинеарностью где N (0, I) нормальное распределение с нулевым средним значением и единичной матрицей ковариаций: E vv T = I.

Упражнение 13.5. Докажите, что в данном примере нормальная сис тема AT AN = AT z относительно МНК-решения xN Rn (здесь n = 2) x имеет вид N +1 xN = (13.21) 11 с общим решением N + 1 xN = x + x = R, N, N, (13.22) 1 где x есть нормальное псевдорешение, единственное МНК -решение с N минимальной нормой среди всех МНК -решений xN, любое число.

Пригоден ли здесь метод нормальных уравнений? Очевидно нет, так как нормальная система (13.21) вырождена.

1. По поводу метода нормальных уравнений. Даже если нормальная сис тема не является вырожденной, она может оказаться плохо обуслов ленной. В данном примере она вырождена: = N. Кроме того, формирование нормальной системы может потребовать хранения и затем одновременной обработки большого массива данных. Напри мер, при обработке данных в астрономии или спутниковой геодезии N (число экспериментальных данных в векторе z) очень велико десятки или сотни тысяч при том, что число оцениваемых параметров n в векторе x сравнительно мало несколько десятков [68]. Единовремен ная обработка сверхбольших массивов данных A z в таких задачах, очевидно, нецелесообразна. Главное, однако, не в этом, а в особенностях численного решения задач с мультиколлинеарностью.

2. По поводу численного решения задачи. Решение нормальной системы в ситуациях, близких к мультиколлинеарности, любым одновременным методом (см. разд. 12, стр. 264) не дает результатов. Так, в данном примере (13.19) метод Хаусхолдера, примененный к нормальным урав нениям (13.21), приводит их к эквивалентному виду N +1 xN =, (13.23) 0 0 с общим решением (13.22), но не позволяет выделить x.

N 13 Устойчивые алгоритмы фильтрации Может ли в подобных задачах помочь последовательная обработка, изло женная выше и принципиально возможная как в информационной форме (см. подразд. 11.6, стр. 244), так и в ковариационной форме (см. под разд. 11.7–11.8 и далее подразд. 13.3–13.8)?

Ai Введем обозначения: A1 = [1 1], z1 = [1], aT = [1 1], Ai+1 =, i aT i zi zi+1 =, i = i, i = 1, 2,.... Тогда AN = A, zN = z. Информационный i+ последовательный алгоритм для i = 1, 2,... имеет вид:

00 i = i1 + ai aT, 0 =, di = di1 + ai i, d0 =.

i 00 N (N + 1) Тогда N =, dN =, и при i = N получаем нормальную систему (13.21) в виде N xN = dN. Отсюда видно, что при таком 0 этот информационный последовательный алгоритм также не решает проблему.

Докажите, что если взять 0 = 2I, то i = Упражнение 13.6.

i + 2 i i(i + 1) = 2, di =, и нормальная система имеет вид i i+ 1 0 i + 1 i(i + 1) i xi = di, следовательно, xi = 1 di = = x.

i i 2 1 2(2i + ) Убедитесь, что прием, примененный в этом упражнении 13.6, фактически совпадает с ридж-оцениванием [18] (гребневое оценивание), т. е. с методом регуляризации Тихонова [67, 81]. Проблема, которая здесь остается, выбор малого ридж-параметра 2. Из этого упражнения видно, что хотя рекуррент ный информационный алгоритм и прост в реализации, и формально дает решение, стремящееся к нормальному псевдорешению x при 0, для яв i ного его получения не избежать решения нормальной системы. При малом 2 это снова создает проблему. Ридж-оценивание ее не решает. Выход можно искать в рекуррентном ковариационном алгоритме МНК.

Упражнение 13.7. Возьмите к качестве первого рекуррентного кова риационного алгоритма МНК скаляризованную форму фильтра Калмана (см. подразд. 13.3) с начальной матрицей ковариации P0 = 1 = 2I при малом (в пределе 0). Получите этот алгоритм в виде:

13.13 Пример задачи с мультиколлинеарностью Таблица 13.1. Скаляризованный фильтр Калмана в задаче регрессионного модели рования с мультиколлинеарностью регрессоров Величина Зависимость от Предел при 0 То же при i = N 2i + 2 i N, i i 2(i 1) + 2 i1 N 1 1 1 1 Ki 1 1 2i + 2 2i 2N i + P0 P0 P i 1 1 1 Pi i i + 2 1 1 1 2i + 2 2 i(i + 1) i+1 N + 1 1 = x xi N 1 1 2(2i + 2 ) 4 Для i = 1, 2,... выполнять:

i = aT Pi1ai + 1, Ki = Pi1 ai /i, (13.24) i Pi = Pi1 Ki aT Pi1, xi = xi1 + Ki(i aT xi1).

i (13.25) i Упражнение 13.8. Докажите результаты работы этого алгоритма, све денные в табл. 13.1.

Таким образом, алгоритм (13.24), (13.25) не требует ни решения систем, ни обращения матриц для задач вида (13.19), (13.20), независимо от раз и при малом 2 он гарантирует выде мерностей и конкретных данных, ление в явном виде решения x в составе обшего МНК-решения (13.22).

N Проблемы выбора малого не возникает;

2 нужно выбирать настолько большим, насколько позволяет диапазон вещественных чисел на компью тере, поскольку параметры алгоритма Ki, Pi обычно быстро убывают и ста билизируются [61].

Упражнение 13.9. Наилучшие результаты дадут в этой задаче не формулы (13.24), (13.25), а их алгебраический эквивалент, получаемый при переходе от Pi к квадратному корню из Pi. Аналогично упражнению 13.8, испытайте другие алгоритмы, а именно: стабилизованный фильтр Калмана– Джозефа из подразд. 13.4, фильтр Поттера из подразд. 13.5, LDLT -фильтр Бирмана из подразд. 13.7, U DU T -фильтр Бирмана из задания на стр. 311, LLT -фильтр Карлсона из подразд. 13.8 или U U T -фильтр Карлсона из зада ния на стр. 311.

13 Устойчивые алгоритмы фильтрации x 2i + i 2i = i+ = 2i 1 ) ) 2( i + + i( = x i+ 0 Рис. 13.1. Сечение поверхности критерия качества J(x) = (z Ax)T (z Ax) (для упражнения 13.10) плоскостью уровня Исследование задачи (13.19) будет неполным, если вы не изучите поверх ность критерия качества J(x) = (z Ax)T (z Ax) и не используете возмож ность преобразования базиса в пространстве Rn, здесь n = 2 размерность оцениваемого вектора x.

Упражнение 13.10. Задайте начальное значение ковариации P0 = = 1 = 2I. Возьмите зависимость Pi от из табл. 13.1. Найдите соб ственные значения матрицы Pi. Убедитесь, что они равны:

1 1 =, 2 =.

2i + 2 Найдите соответствующие собственные векторы. Убедитесь, что они равны:

1 v1 =, v2 =.

Постройте график на плоскости x = (x1, x2 ) и дайте необходимые обоснова ния к его построению (рис. 13.1).

13.14 Задание на лабораторный проект № = Введите преобразование базиса по формуле x x.

2 1 Убедитесь, что при этом P0 = 2I выводится из P0 = 2I. Докажите, что стандартный ковариационный алгоритм дает на i-м шаге значения 2 i(i + 1) 1 1 0, xi = i = i, Ki =, Pi = 2i + 2 0 0 2 2(2i + 2 ) вместо значений во втором столбце табл. 13.1 на стр. 303. Переход к коорди T соответствует методу главных компонент [18]. Видно, натам x что первая компонента выделяется (оценивается) со все возрастающей точ = i(i + 1)/( 2(2i + 2)) (i + 1)/(2 2), ностью (по мере i, 0):

а оценка второй компоненты остается равна нулю: = 0. Проверьте, что ( + )/ обратное преобразование дает правильный результат: x (i + 1)/4 и x2 ( )/ 2 (i + 1)/4.

13.14 Задание на лабораторный проект № Введение В данном лабораторном проекте мы предлагаем к изучению современ ные методы построения численно устойчивых и экономичных по затратам ресурсов ЭВМ алгоритмов метода наименьших квадратов (МHК), решаю щих актуальную задачу последовательного обновления оценок по измере ниям. Как уже отмечалось, эта задача возникает во многих приложениях, например, при оценке состояния (элементов движения) объекта по после довательно поступающим данным наблюдения (см. выше подразд 13.12), при подгонке параметров модели под результаты продолжительных экспери ментов, проводимых в физике, астрономии, экономике, бизнесе, социологии, психологии и в других областях, где выявление регрессий, анализ и прогно зирование тенденций опирается на включение в обработку данных наблюде ния по мере их поступления, для того, чтобы постепенно идентифицировать модель объекта (процесса) или уточнять параметры модели.

1. Задание A. Спроектировать, отладить и продемонстрировать в действии про грамму решения несовместной системы уравнений Ax b, A = A(m, n), 13 Устойчивые алгоритмы фильтрации m n, rank A = n, в смысле наименьших квадратов в соответствии с вашим вариантом последовательного алгоритма (см. ниже подразд. 13.15).

Б. Оценить результаты решения по трем показателям:

(1) погрешность (абсолютная и относительная) решения;

(2) затраты основной памяти компьютера на хранение данных, необходи мых только для заданного алгоритма;

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

Эти показатели определить в зависимости от следующих параметров задачи:

(1) размерность задачи, т. е. число неизвестных n;

(2) степень переопределенности задачи, т. е. число p, указывающее, во сколько раз число уравнений m больше числа неизвестных, m = pn;

(3) степень несовместности системы, т. е. вещественное положительное число c, показывающее среднеквадратическое значение элементов слу чайного разностного вектора d = b A, где x нормальное псевдоре x шение системы Ax b, относительно среднего (единичного) значения;


(4) способ генерации матрицы A.

В. Провести вычислительный эксперимент, используя в нем:

(1) десять значений n, n = 1, 2,..., 10;

(2) три значения p, p = 10, 100 и 1000;

(3) три значения c, c = 1/10, 1 и 10;

(4) четыре способа генерации матрицы A (см. п. 2 в подразд. 13.14).

Результаты эксперимента вывести на экран в виде следующей таблицы:

Вычислительный эксперимент p=(значение), c=(значение), A=(способ) Погрешность Память Число операций n абсолютная относительная КБайт теоретическое реальное При подсчете числа операций учитывать: сложение, умножение, деление и извлечение квадратного корня. К отчету о работе приложить расчетные 13.14 Задание на лабораторный проект № формулы числа операций отдельно по этим видам операций и их сумму.

В таблицу выводить только суммарное число операций.

Г. Выполнить отладку программы и продемонстрировать результаты отладки на решении следующей тестовой задачи [82]:

.

ai1. ai2, Ax b, A = A(m, 2) = i = 1, 2,..., m ;

.

m = 4, 8, 12, 16, 20, 24, 28, 32, 36, 40 ;

ai1 = wi = sin(2i/m), ai2 = wi1 ;

bT = (b1,..., bm), bi = 2 cos(2i/m).

Известно (стр. 267), что решением этой задачи является вектор xT = (1, x2), x x2 = 2 cosec(2/m).

x1 = 2 ctg(2/m), (13.26) Для демонстрации процесса отладки вывести на экран еще одну таблицу:

Отладка программы Погрешность m абсолютная относительная Д. Во всех случаях для оценки абсолютной погрешности использовать норму вектора e = x x вида = max |ei |, e i где x вычисленное решение, а относительную погрешность определять по формуле = e / x, в которой x точное МHК-решение (нормальное псевдорешение задачи).

В отладочной (тестовой) задаче по п. Г решением является вектор (13.26), а в задачах вычислительного эксперимента по п. В вектор x = (1, 2,..., n).

2. Генерация задач для вычислительного эксперимента Как отмечено, для этих задач точное МHК-решение следует задать в виде вектора x = (1, 2,..., n). Затем следует сгенерировать матрицу A (см. ниже) и образовать вектор = A. К нему нужно добавить случайный вектор b x число (см. выше п. 1), а N (0, 1) d = c, где c вектор случайных чисел (от подпрограммы псевдослучайных чисел), взятых из стандартного 13 Устойчивые алгоритмы фильтрации нормального распределения с нулевым средним значением и единичной дис персией. В результате получаем вектор b = A + d для системы уравнений x Ax b.

Для генерации матрицы A необходимо предусмотреть один из четырех способов:

Способ 1 матрица A = A(m, n) заполняется случайными числами из равномерного распределения в диапазоне [100, +100]. Условно запишем это так:

A = [ Random(m n) ].

Способ 2 верхняя часть матрицы A, а именно, ее первые n строк, заполняются по способу 1, а остальные строки образуют подматрицу, в кото рой только первые q столбцов, где q = [n/2] ближайшее снизу целое число к n/2, заполняются как в способе 1, а все остальные столбцы нулевые.

Таким образом, матрица имеет вид:

Random(n n) A= Random Способ 3 первая часть матрицы A должна формироваться как в спо собе 2, а остальная часть образуется располагающимися последовательно вниз блоками из единичных матриц I размера n n:

Random(n n) I(n n) A= ··· I(n n) Способ 4 верхняя часть матрицы A должна формироваться как в способе 2, остальная же часть строится подобно способу 2, но ненулевые q столбцов нижней подматрицы заполняются не случайными числами, а рас полагающимися последовательно вниз блоками из единичных матриц I раз мера (q q):

Random(n n) I(q q) A= ··· I(q q) Замечание 13.18. Hе нужно генерировать всю матрицу A, так же как и весь вектор b и весь вектор d, единовременно. Матрицу A нужно 13.15 Варианты задания на лабораторный проект № генерировать построчно, а векторы b и d поэлементно. Hапример, если a = (a1, a2,..., an ) есть текущая строка матрицы A, то z = aT x + v, где T z текущий элемент вектора b, v текущий элемент вектора d, v = c, N (0, 1) текущее случайное число из стандартного нормального рас пределения. Таким образом, последовательно генерируемые данные (aT, z) нужно вводить в алгоритм решения, а также использовать в нем значе ние r = c2, имеющее смысл дисперсии ошибки измерения вектора = A, b x исключительно последовательно.

13.15 Варианты задания на лабораторный проект № Общее число вариантов составляет 26 (учитывая подварианты разли чия в методах ортогонализации в некоторых из вариантов).

Замечание 13.19. Для всех ковариационных алгоритмов (варианты 1–13) в качестве начальных значений можно взять: x0 = 0, P0 = (1/2)I, 0.

Вариант 1. Стандартный ковариационный алгоритм (Калмана).

Найдите его на стр. 274, подразд. 13.3.

Вариант 2. Стабилизированный ковариационный алгоритм (Джозефа).

Найдите его на стр. 275, подразд. 13.4:

(ii) Обработка наблюдений (очередные данные z = aT x + v):

= aT P a + r, K = P a/, P = (I KaT )P (I aK T ) + rKK T. (13.27) Замечание 13.20. Вычислительные затраты существенно зависят от способа программирования выражений. Hапример, выражение (13.27) для P может быть запрограммировано в следующей последовательности:

W1 = I KaT, (n n)-матрица (n n)-матрица W2 = W1 P, P = W2W1 + r(KK T ) T 13 Устойчивые алгоритмы фильтрации или, эквивалентно, в виде:

v1 = P a, n-вектор P1 = P Kv T, (n n)-матрица v2 = P1 a, n-вектор P = (P1 v2 K T ) + (rK)K T, и в обоих способах можно экономить вычисления, благодаря симметрии матрицы P. Однако второй способ имеет на порядок меньше вычислений:

в первом способе выполняется (1, 5n3 + 3n2 + n) умножений, а во втором только (4n2 + 2n) умножений.

Вариант 3. Квадратно-корневой ковариационный алгоритм Поттера.

Найдите его на стр. 276, подразд. 13.5, в следующем виде.

(i) Инициализация (начальные значения x0, P0 ):

1/ x := x0, S := P0.

(ii) Обработка наблюдений (очередные данные z = aT x + v):

f = S T a, = f T f + r, = 1/(1 + r/), S = S Kf T, K = Sf /, x = x + K(z aT x).

(iii) Экстраполяция (между повторениями этапа (ii)):

S := S, x := x.

Замечание 13.21. Вариант 3, в котором вместо S := S предусмот рена процедура триангуляризации S := triang S, дает следующие версии:

нижние треугольные (S L), или – Версия 3.1 : обе матрицы, S и S, верхние треугольные (S U ).

– Версия 3.2 : обе матрицы, S и S, Именно для этого этап (iii) должен содержать, вместо S := S, процедуру триангуляризации S := triang S, матрицы S. Возможны четыре алгоритма этой процедуры: (1) отражения Хаусхолдера, (2) вращения Гивенса, (3) клас сическая Грама–Шмидта ортогонализация и (4) модифицированная Грама– Шмидта ортогонализация (см. лабораторный проект № 6). Соответственно 13.15 Варианты задания на лабораторный проект № этому, всего имеем 8 подвариантов для указанного варианта 3, сведенных в следующую таблицу:

SL SU triang Хаусхолдер 3.1.1 3.2. Гивенс 3.1.2 3.2. ГШО 3.1.3 3.2. МГШО 3.1.4 3.2. Вариант 4. Факторизованный LDLT ковариационный алгоритм (Бир мана). Найдите его на стр. 279, подразд. 13.7.

Замечание 13.22. Рациональное программирование этого LDLT факторизованного алгоритма должно экономить память компьютера. Здесь можно записывать D поверх D и столбцы L поверх столбцов L (при этом нужна только поддиагональная часть этих столбцов.

Вариант 5. Факторизованный U DU T ковариационный алгоритм (Бир мана). Выведите его самостоятельно аналогично алгоритму варианта 4 [15].

Вариант 6. Факторизованный U U T ковариационный алгоритм Карлсона. Выведите его самостоятельно, опираясь на вывод LLT алго ритма, данный в подразд. 13.8. Получите следующий результат [15]:

(i) Инициализация (начальные значения x0, P0 ):

1/ x0 := x0, U := P0.

(ii) Обработка наблюдений (очередные данные z = aT x + v):

f = UTa, f T = (f1,..., fn ),.

U11f1. 0... T 0 = r, K1 =.

.

Для j = 1,..., n выполнять j = j1 + fj2, j = (j1 /j )1/2, j = fj /(j j ), Ujj = j Ujj, Uij = j Uij j Kj (i), i = 1, 2,..., j 1, j = 1, Kj+1(i) = Kj (i) + fj Uij, i = 1, 2,..., j.

По завершении цикла выполнить x := x + K(z aT x).

K = Kn+1/n, 13 Устойчивые алгоритмы фильтрации Вариант 7. Факторизованный LLT ковариационный алгоритм Карлсона. Найдите его на стр. 285, подразд. 13.8.

Вариант 8. Редуцированный фильтр Бирмана. Найдите его на стр. 287, подразд. 13.9.

Вариант 9. Редуцированный стандартный ковариационный алгоритм Бар-Ицхака–Калмана. Найдите его на стр. 288, подразд. 13.10.

Замечание 13.23. В этом варианте используются только способы и 4 генерации матрицы A (см. выше п. 2 в подразд. 13.14). Инициализация и обработка первых n наблюдений выполняются, как в варианте 1. После этого каждое очередное наблюдение имеет вид z = (a(q) )T x(q) + v, так как T. T. T aT = a(q). 0 и xT = x(q). x(s), где s = n q. Вектор x, так..

же как и соответствующие ему векторы x, x, x, разбит на две части вида:

x(q) размерности q и x(s) размерности s. Запишем алгоритм Бар-Ицхака– Калмана для варианта 9:

T (ii) Обработка наблюдений (очередные данные) z = a(q) x(q) + v:

T = a(q) P (qq) a(q) + r, K (q) = P qq a(q) /, T P (qq) = P (qq) K (q) a(q) P (q), T x(q) = x(q) + K (q) (z a(q) x(q) ), K (sq) = P (sq) (P (qq) )1, P (sq) = K (sq) P (qq), P (ss) = P (ss) K (sq)(P (qq) P (qq) )(K (sq))T, x(s) = x(s) + K (sq) ((q) x(q) ).

x Здесь одинарный верхний индекс в скобках указывает размерность вектора, а двойной верхний индекс в скобках указывает размер матрицы. Матрицы P и P из варианта 1 здесь разбиты на блоки по схеме:

P (qq) P (qs) P = PT.

P=, P (sq) P (ss) Вариант 10. Редуцированный стабилизированный ковариационный алгоритм (Бар-Ицхака–Джозефа). Как и в варианте 9, здесь применяется декомпозиция, т. е. разбиение на блоки векторов и матриц, что выделяет 13.15 Варианты задания на лабораторный проект № редуцированную часть алгоритма вида (13.12) и специфическую часть, напо добие (13.13). Для инициализации и обработки первых n измерений исполь зуется алгоритм варианта 2.

Вариант 11. Редуцированный квадратно-корневой алгоритм Бар Ицхака–Поттера, с нижнетреугольным разложением, S L. За основу берется алгоритм варианта 9, но в нем редуцированная стандартная часть заменяется на редуцированную часть размера (q q). Формулы специфиче ской части заменяются, соответственно, на следующие выражения:

K (sq) = L(sq) (L(qq) )1, L(sq) = K (sq) L(qq), L(ss) = L(ss).

Этап экстраполяции выполняется, как в алгоритме Поттера, вариант 3, при этом L(qq) SL= L(sq) L(ss) и L(qq), L(ss) нижние треугольные матрицы (с верхним знаком или ).

Для инициализации и обработки первых n измерений используется алгоритм варианта 3.

Вариант 12. Редуцированный квадратно-корневой алгоритм (Бар Ицхака–Бирмана–Медана), с P = LDLT -разложением по типу алгоритма Бирмана (см. вариант 4).

Здесь используются декомпозиции вида:

P (qq) (P (sq) )T L(qq) 0 Dq P=, L=, D=, 0 Ds P (sq) P (ss) L(sq) L(ss) где L нижние треугольные с единичной диагональю матрицы, D диа гональные матрицы. Действует замечание 13.23 после Варианта 9, но для инициализации и обработки первых n наблюдений используется алгоритм варианта 4. Для обработки каждого из остальных измерений применяется следующий алгоритм:

(1) Выполнить алгоритм из варианта 4, но в применении к векторам раз мерности q и матрицам размера (q q), т. е. по данным z, a(q), L(qq), D(q), x(q) найти L(qq), D(q), x(q).

13 Устойчивые алгоритмы фильтрации (2) Вычислить K (sq) = L(sq) (L(qq) )1, L(sq) = K (sq) L(qq), L(ss) = L(ss), D(s) = D(s), x(s) = x(s) + K (sq) ((q) x(q) ).

x Вариант 13. Редуцированный квадратно-корневой ковариационный ал горитм (Бар-Ицхака–Карлсона). Действует замечание 13.23 после Вари анта 9, но для инициализации и обработки первых n измерений используется алгоритм варианта 7. Для обработки остальных измерений применяется сле дующий алгоритм:

(1) Выполнить алгоритм из варианта 7 (его следует вывести самостоя тельно, наподобие алгоритма варианта 6), но в применении к векторам размерности q и матрицам размера (q q), т. е. по данным z, a(q), L(qq), x(q) найти L(qq) и x(q).

(2) Вычислить все остальные матрицы и вектор x(s).

Вариант 14. Стандартный информационный алгоритм (см. стр. 245).

(i) Инициализация (начальные значения x0, 0 ) :

..

.d 0. d d0 = 0 x0, :=.

..

(ii) Обработка наблюдений (очередные данные z = aT x + v):

...

.d.d aT. z := +a /r.

...

(iii) Выдача результата: x = 1 d.

В качестве начальных значений рекомендуется взять x0 = 0, 0 = 0.

Вариант 15. Квадратно-корневой информационный алгоритм. См. под разд. 7.3, стр. 110, формулу (7.7), здесь ее рекуррентная версия (13.28).

(i) Инициализация (начальные значения R0, x0) :

z0 = R0 x0;

R0 z0 = R z0.

13.15 Варианты задания на лабораторный проект № (ii) Обработка наблюдений (очередные скалярные данные z = aT x + v):

Rj zj Rj1 zj = Tj, j = 1, 2,..., m, (13.28) aT 0e z где Tj ортогональная матрица, которая выбирается так, чтобы при Rj вести к верхнетреугольному виду расширенную матрицу,j aT номер очередного измерения, все матрицы R здесь имеют размер (nn), при этом все Rj, j 1, верхние треугольные.

(iii) Выдача результата: x = Rj zj.

Hачальные значения рекомендуется взять в виде x0 = 0, R0 = 0. В каче стве преобразований Tj возможны четыре процедуры, указанные в описании варианта 3. Таким образом, всего имеем 4 разновидности данного варианта:

Tj вариант Хаусхолдер 15. Гивенс 15. ГШО 15. МГШО 15. Замечание 13.24. В каждом варианте лабораторного проекта № необходимо выполнить проверку всех положений для задачи с мультикол линеарностью, которые изложены в подразд. 13.13 на стр. 300–305.

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

14.1 Задача оценивания Перепишем постановку задачи оценивания, определенную общими выра жениями (13.1), с указанием дискретного времени не в скобках, а в индексе.

Пусть дискретная динамическая система описывается уравнениями:

xt+1 = t xt + Bt ut + Gt wt, (14.1) zt = Ht xt + vt, (14.2) где xt n-мерный вектор состояния системы, zt доступный m-мерный вектор измерений, ut детерминистский r-мерный вектор (входное управ ляющее воздействие), t = 0, 1,..., дискретное время. Матрицы t Rnn, Bt Rnr, Gt Rnq, Ht Rmn, Qt Rqq, Qt 0 и Rt Rmm, Rt 0 известны, могут быть параметризованы и полностью характери зуют систему (14.1), (14.2). Последовательности {w0, w1,... } и {v1, v2,... } независимые нормально распределенные последовательности шумов с нуле выми средними значениями и ковариационными матрицами Qt и Rt, соот ветственно. Без потери общности считаем, что {w0, w1,... } и {v1, v2,... } не зависят от начального вектора состояния системы x0, распределенного 14.2 Блочные алгоритмы в исторической перспективе по нормальному закону с математическим ожиданием x0 и ковариацией P0, т. е. x0 N (0, P0 ).

x Задача оценивания заключается в получении оценки неизвестного век N тора xt, t = 0, 1,... по доступным наблюдениям Z1 = {zt, t = 1, 2,..., N }, содержащим информацию о векторе xt. Если N = t 1, оценка xt век- тора называется экстраполяционной оценкой, точнее, оценкой одношагового предсказания. Если N = t, оценка xt вектора называется отфильтрованной оценкой. Если N t, оценка xt вектора называется сглаженной оценкой.

Задачи сглаживания в этом пособии не рассматриваются. Они подробно представлены в [61], где можно видеть, что в составе алгоритмов сглажи вания присутствуют оптимальные оценки от фильтра Калмана.

Оптимальные оценки это те, которые наилучшим образом (в зара нее определенном смысле) соответствуют истинному значению вектора xt.

Если критерий качества оценивателя Jt определен условным математиче ским ожиданием как Jt E xT xt Z1 = E (xt xt )T (xt xt ) Z1, где N N t E {·} оператор математического ожидания, xt xt xt погрешность, то оцениватель xt, минимизирующий критерий Jt, называется оптимальным в среднеквадратическом смысле оценивателем.

14.2 Блочные алгоритмы в исторической перспективе Дискретный стандартный ковариационный фильтр (СКФ) Калмана дается уравнениями (13.3). Их вывод широко известен, см. например, [61, 113]. В обозначениях системы (14.1), (14.2) эти уравнения выглядят несколько более обозримо в следующей записи:

+ Этап экстраполяции: t = 0, 1,... ;

P0 = P0, • оценка: x = t x+ + Bt ut, t+1 t (14.3) ковариация ошибки: Pt+1 = t Pt+ T + Gt QtGT. (14.4) t t Этап обработки измерений (фильтрация): t = 1, 2,..., • Kt = Pt HtT (HtPt HtT + Rt )1, (14.5) оценка: x+ = x + Kt (zt Ht x ), t t t (14.6) ковариация: Pt+ = Pt Kt Ht Pt, (14.7) где все данные берутся из модели (14.1), (14.2): t, Gt, Bt, Ht, Qt, Rt известные матрицы-параметры линейной дискретной динамической 14 Ортогонализованные блочные алгоритмы системы, ut вектор управления, zt доступный вектор наблюдений, xt вектор состояния системы с начальным значением x0 N (0, P0 ). x Соединяя в (14.3)–(14.7) два этапа в один, для t = 0, 1,... получим (14.8) так называемую предиктивную форму фильтра Калмана:

и (14.9) оценка: x = t x + Pt HtT (Ht Pt HtT + Rt )1(zt Ht x) t+1 t t (14.8) с {·} = x0 при t = 0 в (14.8), поскольку {·} в (14.8) выражает оценку (14.6), отфильтрованную в результате обработки измерения, но измерений в момент t = 0 еще нет (они начинаются с момента t = 1);

ковариация: Pt+1 = Gt Qt GT + t Pt Pt HtT (Ht Pt HtT + Rt )1Ht Pt T t t (14.9) с {·} = P0 при t = 0 в (14.9), поскольку {·} в (14.9) выражает ковариацию (14.7) оценки (14.6), отфильтрованной в результате обработки измерения, но измерений в момент t = 0 еще нет (они начинаются с момента t = 1).

Уравнение (14.9) есть уравнение Риккати [109] относительно Pt+1.

Замечание 14.1. Видно, что в (14.8) отсутствует слагаемое Bt ut.

Оно опущено в силу предположения ut 0. В любой момент оно может быть добавлено, если ut = 0. Такой прием убрать это слагаемое (для простоты записей) или добавить его (в нужный момент) можно применять всегда.

Квадратно-корневая реализация уравнений (14.3)–(14.9) c использова нием разложений вида Pt± = St±(St± )T рассмотрена выше (разд. 13). Там уже возникала необходимость процедуры триангуляризации, т. е. процедуры приведения матрицы St± к требуемому треугольному виду. Эту идею (при менительно к задаче фильтрации) предложил и разработал Schmidt, S.F.

[129] в 1970 г. В частности, Schmidt показал, что уравнение (14.4) для пред сказания матрицы ковариации ошибки оценивания на этапе экстраполяции может быть заменено на эквивалентное уравнение (St+)T (St+1)T =T, (14.10) t t QT GT где St+1, St+ нижние треугольные матрицы в представлениях Pt+1 = = St+1(St+1)T, Pt+ = St+ (St+)T, Qt = Qt QT. Он также разработал алгоритм t построения ортогонального преобразования T, приводящего к требуемому треугольному виду матрицу, стоящую в правой части формулы (14.10). Эта процедура известна как процедура ортогонализации Грама–Шмидта.

В 1971 г. Kaminski, P.G. предложил новые модификации ковариацион ных и информационных типов квадратно-корневых методов фильтрации, 14.2 Блочные алгоритмы в исторической перспективе обладающие рядом преимуществ перед ранее известными алгоритмами [103].



Pages:     | 1 |   ...   | 4 | 5 || 7 | 8 |
 





 
© 2013 www.libed.ru - «Бесплатная библиотека научно-практических конференций»

Материалы этого сайта размещены для ознакомления, все права принадлежат их авторам.
Если Вы не согласны с тем, что Ваш материал размещён на этом сайте, пожалуйста, напишите нам, мы в течении 1-2 рабочих дней удалим его.