Содержание к диссертации
Введение
Глава 1. Задача параметрической идентификации и современные ор тогонализованные блочные методы 34
1.1. Вопросы параметрической идентификации дискретных линейных стохастических систем 35
1.2. Численная неустойчивость стандартной формы дискретного фильтра Калмана 50
1.3. Классы численно устойчивых реализаций 55
1.4. Квадратно-корневые ортогонализованные блочные алгоритмы 62
1.5. Блочные алгоритмы на основе методов взвешенной ортогонали-зации 76
1.6. Вычислительные аспекты ортогонализованных блочных алгоритмов 92
1.7. Заключение и выводы к Главе 1 103
Глава 2. Алгоритмическое вычисление производных в ортогональных преобразованиях параметризованных матриц 106
2.1. Предпосылки и постановка задачи 107
2.2. Вычисление производных в матричных ортогональных преобразованиях на основе классических методов ортогонализации 112
2.3. Вычисление производных в матричных ортогональных преобразованиях на основе методов взвешенной ортогонализации 140
2.4. Заключение и выводы к Главе 2 155
Глава 3. Адаптивные ортогонализованные блочные методы 162
3.1. Предпосылки и постановка задачи 163
3.2. Общие принципы построения адаптивных ортогонализованных блочных алгоритмов на основе матричных ортогональных преобразований 168
3.3. Адаптивные ортогонализованные блочные алгоритмы 170
3.4. Решение проблемы Дж. Бирмана о вычислении градиента функции правдоподобия в терминах ортогонализованного UD-фильтра 213
3.5. Вычислительные аспекты адаптивных ортогонализованных блочных алгоритмов 225
3.6. Заключение и выводы к Главе 3 237
Глава 4. Развитие метода ВФК для параметрической идентификации с применением адаптивных ортогонализованных блочных фильтров 241
4.1. Предпосылки и постановка задачи 242
4.2. Применение метода ВФК для параметрической идентификации дискретных линейных стохастических систем 247
4.3. Сравнение двух подходов к построению ВФК. Решение проблемы перепараметризации в адаптивном фильтре 259
4.4. Внедрение в метод ВФК адаптивных ортогонализованных блочных алгоритмов 273
4.5. Численные примеры 283
4.6. Заключение и выводы к Главе 4 287
Глава 5. Комплексы программ для исследования и решения задач параметрической идентификации 291
5.1. Предпосылки 292
5.2. «Программный комплекс для параметрической идентификации дискретных моделей стохастических систем PIAAF 1.0» 293
5.3. Параметрическая ВФК-идентификация дискретных моделей суточной динамики температуры человека 297
5.4. «Программа для идентификации параметров в стохастических линейных системах ISLSP v1.1» 318
5.5. Учебный программный комплекс «Стохастические модели, оценивание и управление СМОУ v.2.0» 325
5.6. Заключение и выводы к Главе 5 334
Заключение 338
Список используемых обозначений 345
Список литературы
- Классы численно устойчивых реализаций
- Вычисление производных в матричных ортогональных преобразованиях на основе классических методов ортогонализации
- Решение проблемы Дж. Бирмана о вычислении градиента функции правдоподобия в терминах ортогонализованного UD-фильтра
- Внедрение в метод ВФК адаптивных ортогонализованных блочных алгоритмов
Классы численно устойчивых реализаций
Рассмотрим многомерную дискретную линейную стохастическую систему б, заданную в пространстве состояний уравнениями Xk+i =Fk(S)xk + Bk(0)uk + Gk(0)wk , (1.1) Zk+i =Hk+i(Q)xk+i + Vk+i, к 0 , (1.2) где к - дискретное время, т.е. Xk означает x(tk); Xk Є Шп - подлежащий оцениванию вектор состояния системы и Zk Є Шт - доступный вектор измерений; вектор uk(ERd- детерминированный входной сигнал. Шумы {wk} и {vk} являются независимыми нормально распределёнными последовательностями с нулевыми средними и ковариационными матрицами Qk(@) 0 и i4(@) 0, соответственно. Предполагаем, что средние значения случайных величин известны, поэтому без потери общности считаем их нулевыми. Начальное состояние хо является гауссовым случайным вектором с характеристиками: среднее жо(6) и ковариационная матрица По(6), т. e. хо Л/(жо(6), По(6)). Начальное значение вектора состояния не зависит от {wk} и {vk}.
Предположим, что дискретная система (1.1), (1.2) параметризована по в, т. е. зависит от вектора неизвестных параметров 9єМр, подлежащего оцениванию. Это означает, что дискретная система 6 известна с точностью до определённых параметров, т. е. элементы системных матриц Fk(Q) Є Шпхп, Вк(0) Є M.nxd, Gk(Q) Є M.nxq, Qk(0) Є Шдхд, Hk(Q) Є Штхп и і4(@) Є Штхт могут зависеть от в. Начальные условия жо(6) Є Шп и По (в) Є Шпхп также могут зависеть от неизвестного векторного параметра в.
Согласно Т. Кайлату [46], дискретная система (1.1), (1.2) имеет стандартное представление в виде дискретной модели, которое соответствует физической модели системы. Его можно определить самыми разными способами. Во многих случаях задание вектора состояния определяется конкретной практической задачей. В случае линеаризации нелинейной системы также приходим к стандартной модели. Кроме сделанных здесь предположений относительно независимости шумов {wk} и {vk}, в литературе рассматривают также более общий случай коррелированных шумов (см. подробное изложение, например, в [46, Kailath et al.]), т.е. Е {wkvj} = SkSkj, где Е {} - символ математического ожидания, а 5kj - символ Кронекера.
Если системные матрицы не зависят от дискретного времени к, тогда стохастическая система, заданная дискретной моделью (1.1), (1.2), будет инвариантной во времени (стационарной системой).
Со стандартным представлением стохастической системы (1.1), (1.2) связано множество разных задач: рекуррентного оценивания и адаптивной фильтрации [46, Kailath et al], [64, Огарков], [33, Фомин], параметрической идентификации [16, Льюнг], стохастического управления [131, Maybeck], [17, Саридис], [132, Mosca] и многих других.
Важной задачей, возникающей часто в различных областях науки и техники, является разработка эффективных алгоритмов параметрической идентификации для математической модели (1.1), (1.2) с учётом высказанных выше предположений. Данная задача заключается в нахождении оценок неизвестных параметров в Є W по известным входным U l = [U Q\UJ\ ... гі _1] г и выходным данным наблюдений Zf = [zj\zj\... 4Г в соответствии с выбранным критерием качества идентификации, который определяется некоторым функционалом v7(0; UQ 1) Zi).
При заданном критерии J задача нахождения оценки неизвестного системного параметра в требует решения задачи нелинейного программирования с ограничениями: @ min = argminl7(G; U0 ,Z± ), (1.3) ЄєР(в) где D(0) - область определения параметра в. Наиболее часто для динамических систем вида (1.1), (1.2) оценивание параметров проводят по методу максимального правдоподобия и наименьших квадратов (см. [133, Astrom] и др.). Для метода максимума правдоподобия запишем отрицательную логарифмическую функцию правдоподобия [134, Schweppe]: N 1 N ITMLF(@; UQ , Zl ) = ln(2-7r) У ln(det Re,k) + &k ёкек (1.4) Для метода наименьших квадратов N JMLS(Q;U -\Z ) = -J2{elR lkekV (1.5) k=l где вектор невязки измерений ек и его ковариационная матрица Re при заданных значениях параметра в вычисляются по известным уравнениям дискретного фильтра Калмана [46, Kailath et al].
Наряду с указанными критериями максимального правдоподобия и наименьших квадратов, в диссертационной работе для решения задач параметрической идентификации применяется критерий качества, определяемый вспомогательным функционалом J(S; U \Z ), впервые предложенным проф. И.В. Се-мушиным [68]. Развитие метода вспомогательного функционала качества (ВФК) представлено в Главе 4. Вспомогательный функционал качества также требует применения дискретного фильтра Калмана.
Таким образом, дискретный ФК является основным математическим инструментом для вычисления оценок вектора состояния дискретной линейной стохастической системы.
В начале 1960-х гг. Рудольф Эмиль Калман [135] предложил своё знаменитое решение задачи оптимальной дискретной фильтрации, и это решение в дальнейшем получило название “фильтр Калмана ”. Первоначально ФК был построен для решения задач в области аэронавтики [3, Bierman and Thornton], однако со временем он стал применяться в самых различных областях науки и техники, таких как эконометрика [4, Carraro], молекулярная биология [5, Shamaiah and Vikalo], метеорология [6, Lange], спутниковая геодезия [7, Brockmann], телекоммуникационные сети [8, Murgu], обработка изображений в реальном времени [136, Piovoso and Laplante] и многих других.
Вычисление производных в матричных ортогональных преобразованиях на основе классических методов ортогонализации
Второй класс ортогонализованных блочных алгоритмов построен на основе известных методов взвешенной ортогонализации (см. раздел А.2 Приложения А). Рассмотрим два типа ортогонализованных блочных алгоритмов: UD-реа-лизации и LD-реализации дискретного фильтра Калмана.
Перечислим основные преимущества таких алгоритмов: устойчивость по отношению к ошибкам машинного округления; отсутствие операции извлечения квадратного корня; избавление от операции матричного обращения; компактность и удобство записи ортогонализованной формы; ориентированность на параллельные вычисления.
Отличие данного класса численно эффективных реализаций дискретного ФК заключается в том, что все они основаны на представлении ковариационной матрицы ошибок оценок вектора состояния Pk в виде произведения CDCT, где С - верхняя треугольная U либо нижняя треугольная L матрица с единицами на диагонали, D - диагональная матрица. Для любой квадратной положительно определённой матрицы такое представление можно получить как результат модифицированного разложения Холецкого [57, Bierman], [43, Grewal and Andrews]. Первой UD-реализацией ФК является последовательный алгоритм Дж. Бирмана [57]. Бирман не только доказал вычислительную эффективность UD-фильтра, но и показал, что при соответствующей программной реализации его алгоритм не сложнее, чем стандартный алгоритм Калмана, см. [169, Thorn and Bierman].
Идея построения ортогонализованных блочных алгоритмов на базе модифицированной взвешенной ортогонализации Грама-Шмидта (MWGS – modified weighted Gram-Shmidt orthogonalization) состоит в следующем. Рассмотрим уравнение вида = 1 + 2. Если существует матрица MWGS-преобра-зования такая, что Ат Ст = UB , где – верхняя треугольная матрица c единицами на диагонали, то г -і \ D\ О А А С o 2c .Di О А т т т Рк = А С = A D\A + С D2C = UDU , О D2 С где {U, D} - верхний треугольный и диагональный факторы в UD-разложении матрицы Pk. Другими словами, данный класс методов использует взвешенную ортогонализацию для обновления матриц {U,D} в соответствии с уравнением Риккати, записанным в виде UDUT = ATD\A + С1D C.
Основной вклад в развитие методов реализации UD-фильтра внесли следующие авторы:
1. Дж. Бирман (G.Bierman) [57] предложил использовать UD-факторы ковариационной матрицы ошибки оценки Pk вектора состояния Xk для обновления разностного уравнения Риккати на этапе обработки измерений в фильтре Калмана.
2. КТорнтон (C.Thornton) [150] применила алгоритм взвешенной ортого-нализации Грама-Шмидта для вычисления UD-факторов ковариационной матрицы на этапе экстраполяции в фильтре Калмана.
3. Т. Кайлат (T. Kailath) и соавторы [53] построили новый класс алгоритмов для реализации дискретного ФК - квадратно-корневые ортогонализован-ные блочные алгоритмы (array algorithms). Отличительной чертой таких алгоритмов является удобная форма представления данных для реализации на ЭВМ, ориентированность на параллельные вычисления [52, Jover and Kailath], [170, Hotop]. Кроме того, ортогональные преобразования обладают численной устойчивостью по отношению к ошибкам машинного округления.
4. В [52, Jover and Kailath] авторы предложили ортогонализованный алгоритм этапа обработки измерений в фильтре Калмана, основанный на быстрых вращениях Гивенса. Было показано, что данный алгоритм эквивалентен алгоритму Бирмана для случая скалярных измерений.
5. В книге [43, Grewal and Andrews] авторами сформулирован одностадийный ортогонализованный алгоритм для UD-фильтра, основанный на модифицированной взвешенной ортогонализации Грама-Шмидта. Последовательный алгоритм Бирмана Пусть на этапе фильтрации в алгоритме CKFM (Алгоритм 1) используют разложения Рк = UPkDPkUpk и Рк\к = UPk]kDpk]kUpklk. Предположим, что элементы матрицы Dpk являются положительными, а матрица Rk - диагональная, т. е. Rk = Diag [г (1),..., rk(m)]. Тогда этап фильтрации алгоритма 1 эквивалентен следующему алгоритму: Алгоритм 13 (SBF - Sequential Bierman Filter). 0. НАЧАЛЬНОЕ ПРИСВАИВАНИЕ: U = UPk, D = DPk, x = xk. 1. m-КРАТНОЕ ПОВТОРЕНИЕ ПРОЦЕДУРЫ СКАЛЯРНОГО ОБНОВЛЕНИЯ. Для j = 1,2,... ,m выполнять: Вычислить векторы: J = [/ь Н, . . . , Jn\ = и ti; v = [v\, V2, . . ., vn\ = Vj. Задать начальные значения: а = г + vifi, d\ = d\r/(i\; К = [v\\0 ... 0] . Для і = 2,..., п выполнять: начало а := а + Vifi] 7 := Vа! di := diafry] А := —/«7! uj := uj + ЛІС; ІІГ := ІІГ + fju ; а := а. конец Вычислить векторы v := -f(z - hTx), x := x + Kv с экстраполяцией между повторениями: U = U; D = D; x = x. 2. ЗАВЕРШАЮЩЕЕ ПРИСВАИВАНИЕ: иРщн = U; DPklk = D; xk\k = x. 3. КОНЕЦ. Здесь h - j-й столбец матрицы Щ; z - j-й элемент вектора zk; г - j-й элемент rj{k) диагональной матрицы ковариаций шума измерений Щ, j = 1,..., т - номер скалярного измерения в составе вектора измерений zk в момент времени к.
Доказательство алгоритма можно найти в [57, Biermanl977] на с. 77-80. Замечание 1.9. Алгоритм 13 не содержит операции извлечения квадратного корня, а работа с треугольными матрицами требует меньшего числа арифметических операций по сравнению с обычными. Однако он неудобен для реализации на многопроцессорных вычислительных системах в силу своей последовательной, а не параллельной структуры.
Решение проблемы Дж. Бирмана о вычислении градиента функции правдоподобия в терминах ортогонализованного UD-фильтра
Ортогональные преобразования в силу своих улучшенных вычислительных свойств широко используются при решении различных задач вычислительной линейной алгебры [171, Воеводин], [172, Самарский и Гулин], [173, Райс], [62, Семушин], [152, Голуб и Ван Лоун]. В теории калмановской фильтрации [46, Kailath et al.], [57, Bierman], [43, Grewal and Andrews] ортогональные преобразования применяются для эффективного вычисления решения матричного разностного уравнения Риккати. Под эффективностью подразумеваем прежде всего устойчивость по отношению к ошибкам машинного округления. Многочисленные исследования данной проблемы содержатся в [152, Голуб и Ван Лоун], [57, Bierman], [45, Verhaegen and Van Dooren], [49, Kaminski et al.], [66, Семушин и соавт.] и др.
В свою очередь, потребность в вычислении производных в матричных ортогональных преобразованиях возникает в задачах вычислительной математики, математической физики, теории управления и др. В задаче параметрической идентификации дискретных моделей стохастических систем подобные вопросы требуется решать при построении устойчивых к ошибкам машинного округления алгоритмов вычисления решения уравнения чувствительности Риккати [67, Gupta and Mehra], [133, A strom].
Идея вычисления производных в матричных ортогональных преобразований с приложением к задаче параметрического оценивания впервые была изложена в работе [58, Bierman et al.]. В дальнейшем эта идея получила своё развитие при разработке класса квадратно-корневых методов параметрической идентификации в работах [123–125, Куликова], [107, Семушин и Цыганова], [94, Цыганова], [83, 85, 86, 90, Цыганова и Куликова].
В других предметных областях потребность в дифференцировании матричных ортогональных преобразований возникает в теории возмущений и управления, в дифференциальной геометрии, при решении таких задач, как вычисление экспонент Ляпунова [174, 175, Dieci et al.], задач автоматического дифференцирования [176, Giles], вычисления численного решения матричного дифференциального уравнения Риккати [177, Kunkel and Mehrmann], вычисления производных высокого порядка в задаче планирования эксперимента [178, Walter].
Рассмотрим прямоугольную матрицу А{9), элементы которой зависят от скалярного параметра GR1. Предположим, что область определения V{9) параметра 9 такова, что для любого 9 Є V{9) существует ортогональное преобразование Q(9)A(9) = R(9), где Q{9) - матрица ортогонального преобразования к верхнему треугольному виду прямоугольной матрицы А(9), R{9) - верхняя треугольная матрица. Далее предположим, что элементы матрицы А являются дифференцируемыми функциями по скалярному параметру 9 и матрица производных AQ = \ Ш \ известна. Требуется вычислить значения элементов матри-цы В!в в заданной точке.
Вопросы существования гладких функций, являющихся элементами матриц Q и R, подробно исследованы в [179, Dieci and Eirola]. Далее будем учитывать следующий теоретический результат:
Предложение 2.1. [179, Proposition 2.3] Если элементы матрицы А полного столбцового ранга имеют к-произеодных, тогда существует QR-разложение матрицы А такое, что элементы матриц Q и R также имеют к-производных.
Пример 2.1 наглядно демонстрирует, что при вычислении производных в матричных ортогональных преобразованиях даже очень простые функциональные зависимости элементов параметризованной матрицы () приводят к сложным формулам, для получения которые требуются определённые вычислительные ресурсы специализированных программных систем, таких как Maple, не говоря уже о ручных вычислениях. Следовательно, такой способ вычисления производных никак не подходит для решения задач реального времени, когда процесс решения состоит из множества итераций, при этом на каждой итерации выполняется ортогональное преобразование матриц большого размера, и требуется вычислять значения производных.
Следовательно, актуальной и важной является задача построения таких вычислительных методов, которые позволяли бы найти в заданной точке значения производных элементов матриц, полученных в результате матричного ортогонального преобразования, не прибегая к операции дифференцирования и к громоздким символьным вычислениям. Точность вычисления и скорость работы таких алгоритмов должна быть сравнима с точностью и скоростью вычисления самого матричного ортогонального преобразования.
Цель данной главы — расширить функциональность методов ортогонали-зации и, таким образом, построить и обосновать новые классы методов вычисления производных в матричных ортогональных преобразованиях разных видов.
Для достижения указанной цели были поставлены и решены следующие задачи:
Построить и обосновать новый класс методов вычисления производных в матричных ортогональных преобразованиях вида = и = , где – параметризованная, прямоугольная в общем случае матрица, зависящая от скалярного параметра , – матрица ортогонального преобразования соответствующих размеров, – верхняя правая, – нижняя левая треугольные матрицы. В качестве базовых методов ортогонального преобразования использовать классические методы триангуляризации: метод отражений Хаусхолде-ра, метод вращений Гивенса, модифицированный метод Грама-Шмидта [152, Голуб и Ван Лоун].
Внедрение в метод ВФК адаптивных ортогонализованных блочных алгоритмов
Рассмотрим прямоугольную параметризованную матрицу А{9) и диагональную матрицу Dw{9), где 9 - скалярный вещественный параметр. Найдём решение Задачи 2.2 о построении нового класса вычислительных методов для нахождения в заданной точке 9 = 9 нижней треугольной матрицы производных Ь#І0=0 и диагональной матрицы производных (Dp)fe\0= по известным параметризованным матрицам А{9) и Dw(6), и полученным в результате матричного LD-разложения Ат = ІВТ (где A1DWA = LD U), нижней треугольной матрице L с единицами на диагонали и диагональной матрице Dp.
Предположим, как и в разделе 2.3.1, что область определения параметра 9 такова, что матрица А имеет полный столбцовый ранг и диагональная матрица Dw 0.
В качестве базового метода для вычисления факторов L и Dp используем модифицированный метод взвешенной ортогонализации Грама-Шмидта, см. Приложение А, раздел А.2. Обозначим А = А{9) и Dw = Dw{9).
Смысл LD-факторизации заключается в том, чтобы по заданной паре матриц {Л, Dw} с помощью процедуры модифицированной взвешенной ортогонализации Грама-Шмидта вычислить пару матриц {Z, D } таких, что А = LB и A DWA = LDpL , (2.60) где А є Rrxs, г s и В є Rrxs - матрица MWGS-преобразования, L є Rsxs -нижняя треугольная матрица с единицами на диагонали. Диагональные матрицы Dw є Wxr, Dp є Msxs удовлетворяют условиям BTDWB = Dp и Dw 0 (см. [66, Семушин и соавт, Лемма 3.2, с. 159]).
Пусть элементы матриц {А , Dw} из (2.60) являются дифференцируемыми функциями по скалярному параметру 9. Рассмотрим преобразование (2.60). При известных матрицах производных А в и (Dw) e значения производных элементов матриц {L, Dp} при 9 = 9 можно вычислить как Ц = L(C0 + C2 + Щ) D/ и (Dp) = 2V0 + Т 2 , (2.61) где величины Со, Т о, йо являются, соответственно, строго нижней треугольной, диагональной и строго верхней треугольной частями матричного произведения ВТDWA JJ T, а матрицы Т 2 и &2 являются, соответственно, диа-гоналъной и строго нижней треугольной частями матричного произведения В1[DW) M.
Доказательство. Доказательство этой леммы аналогично доказательству леммы 2.7. Обозначим лишь некоторые особенности. Используем представление А = BU, где U = 17. Кососимметрическую матрицу (Т )ТТ представим в виде (TQ) Т = С[ — С\. Далее используем представления ВтDWAIJJ T = Со + Т о + Ыо , о В1(DW) B = С,2 + Т 2 Все рассуждения аналогичны лемме 2.7. гп
Теоретические результаты, полученные в п. 2.3.1 и 2.3.2, дали возможность построить алгоритмы вычисления производных в матричных MWGS-преобра-зованиях. Сначала рассмотрим UD-разложение параметризованной прямоугольной матрицы А є Wxs, r s, и диагональной матрицы Dw{9) є Wxr, элементы которых являются дифференцируемыми функциями по скалярному параметру 146 9, где 9 Є Т (9). Предположим, что аналитические выражения для элементов матриц производных А в = у - и (Dw) e = -4 - известны. Рис. 2.7 содержит псевдокод алгоритма 2.7, построенного по лемме 2.7. Алгоритм 2.7: Diff.UD Входные данные: А(9) є Rrxs, Dw{9) є Wxr, А в = { } , (В,Ув = {Щ,9 = 9. Результаты: U Є Rsxs, Dp є Rsxs, Щ, {Dp) . 1 начало 3 4 5 6 7 8 9 вычислить А «— А(9)\в=$, А «— А в\в=$,; вычислить Dw - Dw(9)\e=, (Dw) - {Dw) e\e=f найти [U, Dp, В] - MWGS-UD(Л!)„,); // см. рис. А.8 найти X «— ВТDWA U T: о разделить X на три части: [Со + Т о + йо] — X; найти Y «— ВТ(DW) B; разделить Y на три части: [Щ + Т 2+ U i\ — У , получить результат: (Dp) «— 2Т о + Т 2, получить результат: U — U [Щ + UQ + 2 ) Do1. конец
Теперь рассмотрим LD-разложение параметризованной прямоугольной матрицы A EWXS, г s, и диагональной матрицы Dw{9), элементы которых являются дифференцируемыми функциями по скалярному параметру 9, где 9 Є Т {9). Рис. 2.8 содержит псевдокод алгоритма 2.8, построенного по лемме 2.8.
Таким образом, алгоритмы 2.7 и 2.8 позволяют вычислить значения производных элементов матриц, полученных в результате UD- или LD-разложения прямоугольной параметризованной матрицы А{9) в заданной точке 9 = 9 є Ш. Особенность представленных алгоритмов заключается в том, что для вычисления производных факторов {7", Dp} или {Z, Dp} не требуется знания значений производных элементов матрицы MWGS-преобразования В.