Автоматика и телемеханика, № 4, 2011
© 2011 г. М.В. КУЛИКОВА, канд. физ.-мат. наук (Технический университет Лиссабона, Португалия)
МАКСИМАЛЬНО ПРАВДОПОДОБНОЕ ОЦЕНИВАНИЕ ЛИНЕЙНЫХ СТОХАСТИЧЕСКИХ СИСТЕМ В КЛАССЕ ПОСЛЕДОВАТЕЛЬНЫХ КВАДРАТНО-КОРНЕВЫХ ОРТОГОНАЛЬНЫХ МЕТОДОВ ФИЛЬТРАЦИИ1
Рассмотрены теоретические и прикладные аспекты построения численно устойчивых адаптивных методов параметрической идентификации линейных дискретных стохастических систем в пространстве состояний. Неизвестные системные параметры, подлежащие оцениванию, могут входить в любые матрицы, характеризующие систему, и в начальные условия. Впервые предлагается класс градиентных методов, разработанный на основе ортогональных квадратно-корневых реализаций дискретного фильтра Калмана с использованием техники последовательной обработки данных. Показано, что алгоритмы данного типа могут эффективно применяться для решения плохо обусловленных задач параметрической идентификации. Построен тестовый пример. Практическая значимость предлагаемого класса методов иллюстрируется на примере решения одной из задач финансовой математики - идентификации многомерной модели стохастической волатильности.
1. Введение
Рассматривается задача параметрической идентификации линейных дискретных стохастических систем методом максимума правдоподобия. С целью нахождения экстремума функционала правдоподобия, как правило, используются градиентные методы оптимизации или методы ньютоновского типа (см., например, [1]). Вычисление функции правдоподобия для исследуемых систем требует применения дискретного фильтра Калмана, стандартная реализация которого не устойчива по отношению к ошибкам округления (см. [1- 5] и др.). В свою очередь, вычисление градиента функционала заключается в дифференцировании уравнений фильтра, что приводит к системе векторных уравнений, известных как уравнения чувствительности фильтра (filter sensitivity equations), а также к системе матричных уравнений чувствительности вида Риккати (Riccati-type matrix sensitivity equations). Решение уравнения Риккати и, как следствие, уравнений чувствительности фильтра этого типа является основной проблемой как с точки зрения вычислительных затрат, так и с точки зрения устойчивости алгоритма по отношению к ошибкам округления [6].
Из теории фильтрации хорошо известно, что квадратно-корневые реализации дискретного фильтра Калмана позволяют обновлять уравнение Риккати в терминах треугольных матриц, полученных путем разложения Холецкого матриц кова-риации ошибки оценивания Рд|д и предсказанияРд+1|д, и тем самым являются более устойчивыми по отношению к ошибкам округления, т.е. позволяют получить большее количество точных цифр после десятичной запятой, чем стандартная реализация фильтра. Заметим, что уже с момента открытия техники калмановской фильтрации развитие численно устойчивых реализаций является одним из важнейших направлений в данной области исследований [7]. Среди достижений уже был
1 Работа выполнена при финансовой поддержке португальского фонда науки и технологии (Fundagao para a Ciencia e a Tecnologia) и софинансировании европейского фонда FEDER (грант No. SFRH/BPD/64397/2009).
4* 99
отмечен класс квадратно-корневых алгоритмов (square-root filtering algorithms), в развитие которых основной вклад внесли работы [8-10]. Следует также выделить UD фильтр Бирмана, предложенный в [11]; быстрые алгоритмы фильтрации (fast Chandrasekhar-type algorithms) [12, 13]; и, наконец, квадратно-корневые ортогональные реализации (array square-root filtering algorithms), сформулированные сравнительно недавно в [14]. Использование численно устойчивых ортогональных преобразований на каждом шаге рекурсии делает последний класс методов наиболее предпочтительным к внедрению и использованию на практике [15]. Более того, алгоритмы этого типа ориентированы на параллельные вычисления (см., например, [16]).
Существование численно устойчивых реализаций фильтра дает надежду на их эффективное использование и в решении тесно связанной задачи параметрической идентификации. Возникает вопрос о существовании класса адаптивных фильтров, которые позволяют одновременно провести идентификацию системы (здесь рассматривается метод максимума правдоподобия), и оценить ненаблюдаемый вектор состояния динамической системы, т.е. осуществить фильтрацию. Другими словами, возможно ли обновлять в терминах квадратно-корневых фильтров не только уравнение Риккати, но и матричные уравнения чувствительности вида Риккати? И если такие методы существуют, то унаследуют ли они преимущества квадратно-корневых алгоритмов? Прежде всего интересен вопрос устойчивости к ошибкам округления.
Целью статьи является попытка ответить на поставленные выше вопросы и впервые построить класс градиентных методов на основе последовательной реализации дискретного фильтра Калмана из [17]. Похожая задача, но с применением квадратно-корневого информационного фильтра была решена в [18], однако устойчивость алгоритма (к ошибкам округления) не была исследована. Позднее, в [19, 20] были построены реализации, основанные на обработке вектора измерений целиком. Однако, для ряда практических задач (медико-биологических, эконометрических, финансовых и др.), когда модели задаются многомерными линейными дискретными стохастическими системами, важно обрабатывать поступающие данные покомпонентно. Решение поставленной задачи осуществляется в пространстве состояний. Как и ранее, неизвестные системные параметры, подлежащие оцениванию, могут входить в любые матрицы-параметры системы и начальные условия в любых комбинациях.
2. Постановка задачи
Рассмотрим линейную дискретную стохастическую систему:
(1) Хк = Fk (в)хк-1 + Gk(9)wk,
(2) Zk = Нк(в)хк + vk,
где Fk(в) е K"x", Gk(e) е RnX9, Нк(в) е Rmxn, Qk(6) е Wxq и Rk(0) е Rmxm, Хк е К" - вектор состояния системы, zk е Кт - доступный для наблюдения вектор измерений. Последовательности [wo,w\,... } и [v\,v2,... } - независимые нормально распределенные белые последовательности шумов с нулевыми средними и ковариационными матрицами Qk(e) > 0 и Rk(e) > 0 соответственно. Кроме того, wk - N(0,Qk(e)) и vk - N(0,Rk(в)) не зависят от х0 - Я(х0, П0(в)), в е К - вектор системных параметров. Для удобства дальнейшего изложения опустим аргумент в при написании Fk, Gk, Нк, Qk, Rk и По, не забывая при этом об их зависимости, и введем следующие обозначения. Будем рассматривать разложение Холецкого вида: A = AT/2A1/2, где А1!2 - верхняя треугольная матрица.
А-1!2 = (A1/2)-1, A-T/2 = (A-1/2)T и (a(i))T - г-я строка матрицы А.
Поскольку матрицы, характеризующие систему (1), (2), неизвестны, то прежде чем фильтр Калмана может быть использован для нахождения оценки ненаблюдае-
мого вектора состояния жд, необходимо решить задачу идентификации системы, т.е. оценить в по доступным результатам наблюдений ^^ = [21,..., 2^]Т. При решении задачи методом максимума правдоподобия необходимо максимизировать функцию правдоподобия, логарифм которой задается формулой из [21]:
1 м
(3) Се (г?) = - ]Г {-? 1п(2тг) - 1п [с1е1(Де,й)] - е£Д;>й} , £=1
где ед = 2д — ДдЖдд— - невязка измерений дискретного фильтра Калмана с матрицей ковариации ДеД = Дд + РД|Д-1 ДТ. Вектора Ждд— и Ждд - предсказанная и отфильтрованная оценки вектора состояния жд системы (1), (2) соответственно. Рдд— и Рдд - матрицы ковариации ошибки предсказания и оценивания, т.е. Рй|й_1 = Е [(жй — Ждд— — Ждд—)Т] и Рй|й = Е [(жй — ждд.)(жй — Ждд.)Т].
Для максимизации Сд , как правило, используются градиентные методы оптимизации или методы ньютоновского типа, которые, в свою очередь, требуют вычисления градиента функционала. Из (3) получаем
дСв = 1 " ( д (1п Ие1(Д.е,й)]) Э | ,=
Вычисление (3), (4) требует применения дискретного фильтра Калмана, стандартная реализация которого задается следующими рекуррентными уравнениями:
(5) Жд.+1|д. = ДдЖдд— + ДдД-Д (2Й — Дйжй|й_1), Жо = Жо,
где Дд = Д%Р^|й_1ДТ и матрица ковариации РД|Д-1 обновляется рекуррентно (разностное уравнение Риккати):
(6) Рй+1|й = рк|к_1^Т + дйсТ — ад-1дТ, -о = По, По > о.
Хорошо известно, что при реализации алгоритма (5), (6) с помощью ЭВМ возникают проблемы устойчивости вычислений по отношению к ошибкам округления. Одна из них - потеря специального вида матрицами Рд^-ь Рдд (положительно определенные, симметричные). Как следствие, были разработаны квадратно-корневые реализации дискретного фильтра Калмана, основанные на разложении Холецко-го. В данном классе методов разложение осуществляется для начальных условий2,
т.е. По = ПТ/2По/2 и далее обновление (6) происходит в терминах треугольных
1/2
матриц РД|Д-1. На последнем этапе фильтрации, т.е. после обработки последнего
измерения вновь имеем Р^+1|№ = Р_Т+21|№Р^/+1|№ и/или Р№|№ = рТ/2Р^. Данная организация вычислений не свободна от влияния ошибок округления, однако она гарантирует положительные элементы на диагонали ковариационных матриц и их симметричность. Кроме того, было показано, что квадратно-корневые методы позволяют сохранить большее количество точных цифр после десятичной запятой, чем стандартная реализация фильтра Калмана, работающая с полными матрицами Рд | Д-1, Рд | д.
На современном этапе исследований предпочтение отдается ортогональным реализациям из класса квадратно-корневых фильтров. Суть этих методов заключается
2 Для начальных условий По разложение Холецкого существует, так как По — вещественная симметричная матрица с положительными ведущими минорами, поскольку По > 0, т.е. строго положительно определенная.
в формировании блочной матрицы, содержащей в себе всю необходимую для данного этапа фильтрации информацию, и последующем приведении этой матрицы к требуемому нижнему или верхнему треугольному виду с помощью ортогональных преобразований (см., например, [14]). Рассмотрим одну из таких реализаций, развитую в [17]. Принцип ее построения заключается в использовании процедуры последовательной обработки данных. Подобный класс методов выгодно отличается тем, что не требует на каждом шаге итерационного процесса обращения матриц. Тот же принцип рекуррентного обращения матриц используется в алгоритмах дл
Для дальнейшего прочтения статьи необходимо приобрести полный текст. Статьи высылаются в формате PDF на указанную при оплате почту. Время доставки составляет менее 10 минут. Стоимость одной статьи — 150 рублей.