CC BY
30
УДК 681.3:529.7 Серышева Ирина Анатольевна,
старший преподаватель кафедры автоматизированных систем ИрГТУ, e-mail: sia_cyber@mail.ru
Хрусталёв Юрий Петрович, к. т. н, доцент кафедры вычислительной техники ИрГТУ, e-mail: khrustalev@istu.irk.ru
АДАПТАЦИЯ МНОГОМЕРНЫХ ДИНАМИЧЕСКИХ СТОХАСТИЧЕСКИХ МОДЕЛЕЙ
I. A. Serisheva, Yu. P. Khrustalev
ADAPTATION OF MULTIDIMENSIONAL DYNAMIC STOCHASTIC MODELS
Аннотация. В работе рассмотрены вопросы адаптации многомерных динамических стохастических моделей, применяемых в задачах оценивания состояния групповых эталонов (эталонов времени и частоты). Особенность объекта состоит в том, что в результате измерений, выполняемых в процессе его функционирования, находятся разности относительных отклонений частоты опорного генератора со всеми остальными, входящими в состав группового эталона. Это не позволяет производить подстройку параметров прогнозирующих моделей, непосредственно сравнивая прогнозы вектора состояния с их «истинными» значениями. В работе рассматривается метод адаптации, основанный на использовании стохастического квазиградиента. Устойчивость адаптивной процедуры обеспечивается применением операции проектирования вектора параметров в область допустимых значений. Предложен алгоритм, позволяющий решить поставленную задачу. Приведены результаты численных экспериментов, иллюстрирующих возможности разработанного метода.
Ключевые слова: прогнозирующие модели; процессы авторегрессии; адаптация; временные ряды; оптимизация; недоопределенные системы.
Abstract. The questions of adaptation of the multidimensional dynamic stochastic models used in the state estimation problems of group standards (time and frequency standards) are considered. Predictive models are constructed on the basis of a priori information about the time series, available to researchers. When obtaining additional data it is expedient to make tuning ofparameters of models. In this paper the method of adaptation based on the use of stochastic quasigradient is considered. Operation of introduction of vector's parameters in area of admissible values allows providing stability of adaptive procedure. The algorithm, allowing to solve the problem is offered. Results of the numerical experiments illustrating possibilities of the developed method are given.
Keywords: predictive models; autoregression processes; adaptation; time series; optimization; underdetermined systems.
Введение
Задача обеспечения единства измерения времени и частоты на территории Российской Федерации возложена на Государственную службу времени и частоты РФ (ГСВЧ). В составе ГСВЧ функционируют Государственный эталон и вторичные эталоны, являющиеся групповыми, что повышает их точность и надежность. Повышение точностных характеристик эталонов связано как с совершенствованием их аппаратурной базы, так и с применением современных методов обработки данных.
1. Методы стохастической аппроксимации в задачах адаптации прогнозирующих моделей
Задачу обработки данных, получаемых в результате функционирования групповых эталонов физических величин, можно рассматривать как задачу оценивания состояния линейных динамических систем с неполной матрицей наблюдений. Примером таких систем может являться групповой эталон времени и частоты, в котором измерения выполняются через определенные интервалы
времени (как правило, этот интервал равен одним суткам). В процессе измерений получаются разности относительных отклонений частоты опорного элемента - водородного стандарта частоты -и всех остальных элементов группового эталона. Понимая под вектором состояния эталона на шаге 5 относительные отклонения частоты водородных стандартов от приписанных им значений, задачу обработки измерительной информации, получаемой в процессе эксплуатации эталона, можно сформулировать в рамках теории оценивания состояния динамических систем [1]. Решение поставленной задачи известно как «фильтр Калма-на». Произведя декомпозицию фильтра, можно получить формулы для оценки состояния опорного элемента у и всех остальных составляющих
группового эталона у [2]
У =
n
Eft (z + У (!)),
(1)
где Si - вес 7-го измерения;
1
Современные технологии. Математика. Механика и машиностроение
п - число элементов в групповом эталоне;
■' 1 - прогноз значения частоты /-го генератора,
вычисляемый на предыдущем такте обработки данных.
В выражении (1) веса .у; полагаются обратно пропорциональными дисперсиям прогнозов. Для простоты в данной работе будем полагать веса всех измерений одинаковыми, равными —.
Прогнозы У:(Д) вычисляются на основе использования прогнозирующих моделей - динамических стохастических моделей, т. е. моделей авторегрессии скользящего среднего (АРСС).
Для согласования индексов в уравнение (1) добавлено «фиктивное измерение»
г = - = -'. В этом случае индексирование
ведется от 1 до п.
Методика построения моделей АРСС по экспериментальным данным широко известна (методика Бокса - Дженкинса) и применяется при решении различных практических задач. В нашем случае использование данной методики осложняется тем, что в процессе эксплуатации эталона имеют место косвенные измерения величин
т. е. отсутствуют исходные временные ряды. Поэтому был разработан метод оценивания параметров моделей АРСС, основанный на минимизации суммы квадратов отклонения результатов измерений от их прогнозов й^(1) [3]. Поскольку прогнозы измерений зависят от параметров прогнозирующих моделей, минимуму функционала
где^;('1) - прогноз 1-го измерения на такте -5, - ■ - длина временного ряда, соответствуют оптимальные оценки параметров моделей АРСС. (В дальнейшем в работе будет рассматриваться обработка данных на одном из тактов измерений, поэтому индекс 5 будет опущен).
При обработке данных в статическом режиме, т. е. в режиме их накопления, использование прогнозирующих моделей позволяет повысить точность оценивания (т. е. снизить алгоритмическую погрешность) примерно на 20-30 процентов [3]. Традиционные методы повышения точностных характеристик эталонов времени и частоты основываются на совершенствовании их аппаратурной базы. Учитывая высокую стоимость аппаратуры, входящей в состав современных эталонов времени и частоты, можно говорить о высокой эффективности предлагаемого в [3] метода.
В режиме динамической обработки данных, т. е. в режиме их обработки в темпе поступления, возникает необходимость решения следующих задач:
1. В случае стационарных процессов при увеличении объема выборки, т. е. при выполнении текущих измерений, желательно использовать получаемую при этом дополнительную информацию для уточнения параметров прогнозирующих моделей.
2. При разладке моделей, когда прогнозирующая модель становится неадекватной реальному процессу, необходимо включить процедуру перенастройки модели.
3. При изменении остаточной дисперсии, характеризующей средний квадрат ошибки прогнозов, требуется соответствующим образом корректировать вес в формуле (1).
Очевидно, что для решения указанных проблем необходимо использовать адаптивные процедуры.
Теоретической основой адаптивных процедур являются методы стохастической аппроксимации, или, в более общей постановке, методы стохастического программирования. Применение метода стохастического квазиградиента для адаптации моделей АРСС рассматривалось в работе [4]. При этом адаптация одномерной модели сводится к построению последовательности точек
Р*-1 = Кр (Р* - РзУз «5 РД (3 )
где ¡Тр - операция проектирования в область допустимых значений параметров модели; р^- - вектор параметров модели;
"- номер шага; .'.-- коэффициент адаптации; - нормирующий множитель; обобщённый градиент.
Процесс (3) сходится при следующем выборе коэффициента рв [4]
В работе [4] рассматривался упрощенный подход к процедуре адаптации. Во-первых, использовалась одномерная модель АРСС, тогда как на самом деле в групповой эталон входит несколько водородных стандартов, т. е. необходимо рассматривать модели многомерных временных рядов. Во-вторых, что значительно важнее для использования предлагаемого метода при работе с реальным объектом, в распоряжении исследователя нет «истинных» значений временного ряда ув,
а есть только результаты измерений, т. е. вектор
ИРКУТСКИМ государственный университет путей сообщения
(в дальнейшем будем опускать индекс понимая, что речь идет о текущем такте обработки данных).
Все вышесказанное приводит к необходимости обобщения результатов, полученных в [4]. С учетом сделанных замечаний сформулируем задачу адаптации моделей временных рядов, используемых в системах обработки данных, получаемых в групповых эталонах.
Итак, на текущем такте имеем: результаты измерений z — у^ — у,;
вектор прогнозов измерении, вычисленный на предыдущем такте г(1)=Х-р, где X" - транспонированный вектор состояния процесса АРСС [6].
где V- оценка величины у на шаге 5 — /: рп - порядок процесса авторегрессии (АР) 1-го элемента группового эталона; цп - порядок скользящего среднего (СС) 1-го элемента;
= : :: : ■ - расширенный вектор параметров моделей АРСС.
Задача адаптации сводится к построению вскторовр5. минимизирующих целевую функцию
F(P)=^=1(_zs-zs(l-)У, (5)
где ij.fl) - функция расширенного вектора параметров р^, т. е. ¿д("1) = т. е. прогноз ¡'-го измерения, вычисленный на предыдущем такте.
Упрощенная структурная схема процесса адаптации приведена на рис. 1.
В развернутом виде выражение (5) будет представлено следующим образом:
_шшш
- ¿-я составляющая вектора состояния процесса АРСС;
Р*; = [фв1> - - вектор пара-
метров модели! -го генератора.
В (5) полагаем ^1=0, если опорным является 1-й генератор. В этом случае в целевую функцию (5) для опорного генератора включается лишь прогноз его частоты.
Тогда, введя в формулу (3) нормирующий коэффициент у± = сведем задачу минимизации
целевой функции 7"Чр) к задаче поиска ее минимального среднего значения [4].
В выражении (3) обобщенный градиент сг- вычисляется по правилу векторного дифферен-¿рк р:
цирования (х3 —
чр
Последовательность точек р5 строится в многомерном пространстве параметров прогнозирующих моделей. При минимизации целевой функции F (выражение (5)) следует иметь в виду, что от параметров моделей зависят лишь прогнозы вектора измерений ^¿(1), вычисленные на предыдущем такте. В общем случае прогнозы на один шаг вперед находятся как разность прогнозов значений частоты опорного и г-го генераторов
Обозначив параметры авторегрессии для го генератора как фц,фР; и, соответственно, параметры скользящего среднего как бц,— получим выражение для 1-й составляющей вектора
Рис. 1. Схема адаптации многомерных прогнозирующих моделей
Современные технологии. Математика. Механика и машиностроение
дР дГ
дР
обобщенного градиента « дГт _ дГ дГ
(р - порядок АР в модели г -го генератора, ц - порядок СС).
Учитывая, что члены, содержащие коэффициенты скользящего среднего, входят в прогнозы со знаком «минус», а члены, соответствующие коэффициентам авторегрессии, - со знаком «плюс», получим выражение для коррекции составляющих вектора параметров авторегрессии прогнозирующей модели для ! -го генератора
где 7Тф - оператор проектирования в область допустимых значений ф;
Ф}311 I = 1. 2, - /-я компонента подвектора коэффициентов авторегрессии 1-й прогнозирующей модели на шаге
У (я-/) г - значение частоты ¡'-го генератора на шаге
= —¿1.1); - ошибка прогноза 1-й составляющей вектора г.
Поскольку истинные значения частоты водородных стандартов не известны, следует в выражении (6) заменить их оценками у\, которые
находятся в соответствии с алгоритмом субоптимальной фильтрации [2]. Коррекция параметров скользящего среднего производится аналогично. При этом учитывается знак, с которым входят в прогноз эти коэффициенты
= + 2 Аг&^-ы), (7)
где т^е - оператор проектирования в область обратимости коэффициентов СС;
~ разность между оценкой величины 'У{3-/;ц
и ее прогнозом, вычисленным на предыдущем такте.
Формулы (6), (7) позволяют вычислять параметры ф]В1, прогнозирующих моделей для
всех составляющих вектора р на такте 5 обработки данных. Для реализации адаптивной процедуры необходимо еще определить коэффициент адаптации рз как функцию такта
Выбор коэффициента адаптации должен проводиться, исходя из требований сходимости процедуры (6)-(7). Если речь идет о сходимости процесса из любой начальной точки пространства параметров, должно выполняться требование [4]
В противном случае можно выбрать р в соответствии с формулой (4).
Уидроу в работе [7] вводит понятие рабочей функции для обозначения функционалов вида (2), а задачу адаптации рассматривает как задачу поиска минимального значения этой функции в случае стационарных процессов. При нестационарности процессов гЕ (а, следовательно, и ;у5) задача
адаптации может рассматриваться как задача «слежения» за точкой минимума характеристической функции, блуждающей в пространстве параметров.
Рассматриваемая задача несколько отличается от задачи настройки антенной решетки, рассматриваемой в [7], так как мы не имеем возможности бесконечно подстраивать прогнозирующие модели. Более того, с целью уменьшения «шума градиента» (термин, введенный Уидроу), когда подстройка коэффициентов продолжается, даже если найдена точка минимума, мы не будем фиксировать величину р3 на некотором заранее выбранном малом значении. Однако в этом случае необходимо решать проблему контроля адекватности моделей с тем, чтобы после обнаружения разладки модели устанавливать коэффициенты Р;
в максимальное значение, тем самым переходя в режим перенастройки параметров.
При моделировании процесса адаптации рассматривалась также линейно убывающая до некоторого малого значения функция Оче-
видного преимущества перед обратно пропорциональной функцией р3 — ^/г такой выбор не продемонстрировал, но при этом возникает дополнительная проблема: выбор коэффициента Ьопре-деляющего скорость убывания функции. Поэтому в дальнейшем коэффициент адаптации выбирался в соответствии с формулой (4).
Как указывалось выше, нормируемый множитель Уз может включать в себя коэффициенты 2
и Однако их влияние можно компенсировать
соответствующим выбором начального значения р о и не учитывать при дальнейшем рассмотрении процесса адаптации. Необходимость ввода в процедуру подстройки вектора параметров р
нормирующего множителя у5 диктуется требованиями перевода разностей векторов (измеренного и прогнозируемого) в диапазон изменения вектора р. Например, при моделировании рассматривались временные ряды, порождаемые процессами авторегрессии первого порядка со значе-
ИРКУТСКИМ государственный университет путей сообщения
ниями коэффициентов фц. лежащими в диапазоне [-1, 1]. Временные ряды у, (а значит, и^) изменялись в гораздо более узком диапазоне. Чтобы преодолеть это несоответствие, в качестве нормирующего множителя на шаге s выбиралась величина, обратная сумме ¿-норм вектора состояния процесса, т. е.
(под i-нормой вектора в вычислительной математике понимается сумма модулей составляющих вектора [8]). В формулу (3), определяющую процесс адаптации на текущем шаге 2. входит оператор тгр, обозначающий операцию проектирования вектора параметров р в области допустимых значений (т. е. параметров АР - в область устойчивости, параметров СС - в область обратимости).
В [9] эта задача решалась, исходя из геометрических представлений в 3-мерном пространстве. Более общее решение может быть найдено методами математического программирования. В случае линейных ограничений, а мы имеем дело именно с такими ограничениями, можно для решения задачи проектирования использовать метод Розена [10]. Рассмотрим подробнее применение этого метода для решения задачи проектирования параметров АР в область устойчивости (т. к. область обратимости параметров СС описывается аналогичными уравнениями, то мы его отдельно не будем рассматривать).
Проекция точки X на множество М есть точка X <= М, ближайшая к X. Под проекцией вектора f понимается вектор Р € М, менее других отличающийся от f. Проективные методы на каждой итерации вовлекают в вычислительный процесс по возможности минимальное количество ограничений в виде равенств (активные ограничения) [10].
Алгоритм проектирования начинает работу в некоторой допустимой точке Х1-*-'. Затем определяется допустимое направление d = — P7f( X), где Р- матрица проектирования. В этом направлении делается шаг длиной Лк, минимизирующий f(X), но в то же время сохраняющий допустимым вектор Х'-'"--1-' = X'"-- + Ad. Если d=0, имеет место либо оптимальная точка, либо некоторое ограничение несущественно.
Математическая формулировка задачи выглядит следующим образом: min/(X)
при АХ < Ь,
ЕХ = с (8)
Матрица А разбивается на две подматрицы Ат = (А^Ар, а вектор Ь - на два подвекто-ра Ьг = (Ь£, Ьр. При этом А]Х = А:Х < Ь2.
Тогда матрица проектирования имеет вид Р = 1 - Мт ■ (МВТ)-1 - м, где М = (А^, ЕП
В случае с уравнениями авторегрессии второго порядка, ограничения записываются следующим образом:
Пусть координаты начальной точки Ф равны (3, 1). Необходимо найти кратчайшее расстояние от этой точки до X, X Ей. Целевая функция имеет вид: тт/(Х) — Оч- ФаЗ^+О^г- Ф^'при
ограничениях АХ < Ь.
Ограничения-равенства в задаче отсутствуют. Ограничения (8) можно записать в виде:
+ < 1
Тогда матрица А и вектор b равны:
При Ф = (3,1) градиент равен
В этом случае активно второе ограничение: При этом
Так как ограничения-равенства отсутствуют, то м = А И d = -Р ■ УГСО = | .
Двигаясь в указанном направлении и выбирая шаг "/. из условия минимума функции Г(Х'-1; + Яd1), получаем ЛП[ГГ=0,2 и, соответственно, новые координаты Хл = [0,1].
Выполняя следующие шаги (шаг 2 и шаг 3), получаем решение задачи Х- =
Нетрудно убедиться, что в этой точке выполняются условия Куна - Таккера, т. е. найдена проекция исходной точки Ф = [3,1].
Современные технологии. Математика. Механика и машиностроение
Рассмотренные выше процедуры позволяют реализовать процесс адаптации прогнозирующих моделей, используемых в задачах оценивания состояния группового эталона по результатам измерений. Приведенный ниже пример иллюстрирует изложенные выше положения.
При моделировании генерировались 3 временных ряда, соответствующие моделям авторегрессии 1-го порядка с коэффициентами
Среднее квадратическое отклонение белого шума, возбуждающего систему, полагалось равным о = 0.1 для всех трех рядов. Ряды формировались в соответствии с ранее апробированной методикой [5], т. е. начальные значения всех рядов полагались равными нулю (V; = 0.,1 = 1,2,3).
К получаемым на следующем такте прогнозам добавлялись члены О; - элементы выборки из белого
шума и т. д. Мы не приводим фрагментов сгенерированных рядов, т. к. они соответствуют типичным процессам авторегрессии. Уже через несколько тактов получаем временные ряды с установившимися характеристиками. Мы ограничились длиной рядов в 250 точек. Полученные выборки экспортировались в ППП 8ТЛТКТ1СЛ, с помощью которого находились оценки параметров моделей, соответствующих сгенерированным рядам. По этим выборкам получены такие оценки «истинных» параметров моделей: Фи^ = 0.827,
0125С, = 0.323 и = 0.527. Затем, с учетом
того что в распоряжении исследователей имелась выборка ограниченной длины (для определенности полагалось, что выборки состоят из 50 членов каждого ряда), вычислялись оценки параметров моделей АР(1): Фн^ = 0-5. Флз.^ = 01 и
С момента к= 50 начинался процесс адаптации, т. е. подстройки коэффициентов : = 1 _ :■), соответствующий формуле (6).
При этом начальное значение шага Ро полагалось 1
т. е. равным единице. Нормирую-
равным т^,
щие множители Уь на каждом шаге принимались равными сумме ¿-норм значений Уы. Оператор проектирования тг^ сводится к ограничению величин фц < 1). Поправочный член для каждой составляющей вектора параметров (т. е. в нашем случае для коэффициентов фц, Ф12 ■ -_:), определялся взвешенной разностью «изме-
ренной» соответствующей составляющей вектора zki и ее прогнозом вычисленным на преды-
дущем такте обработки данных. Таким образом, мы решаем задачу пошаговой минимизации суммы квадратов ошибок прогнозов вектора измерений, являющейся функцией параметров используемых моделей.
На рис. 2 и 3 приведены графики изменения параметров авторегрессии для каждого из трех рядов и средней суммы квадратов ошибок от прогнозов.
Рис. 2. Процесс адаптации коэффициентов авторегрессии
о.::
555 0.005
1
V /V А,
50
60
30
Рис. 3. Усредненная сумма квадратов остатков прогнозов вектора измерений от «измеренных» значений
В процессе выполненных исследований была проделана серия численных экспериментов. Вполне осознавая, что этого явно недостаточно для того, чтобы говорить о всеобъемлющем статистическом моделировании процесса адаптации, все же можно сделать ряд важных выводов, позволяющих определить направления дальнейших исследований.
1. Как видно из рис. 2, процесс адаптации, двигаясь в правильном направлении, все же не позволяет достичь истинных значений коэффициентов АР, т. е. значений, получаемых по полной выборке при минимизации соответствующего функционала и с использованием методов матема-
ИРКУТСКИМ государственный университет путей сообщения
тического программирования [6]. При этом задача минимизации целевой функции (2) успешно решается (рис. 3). Данное противоречие только кажущееся. В фундаментальной работе [11] приведены графики линий равного уровня для сумм квадратов остатков, получающихся при подгонке коэффициентов модели одного временного ряда с двумя параметрами. Из этих графиков видно, что задачи идентификации и «подгонки» коэффициентов авторегрессии скользящего среднего методом наименьших квадратов имеют множество почти эквивалентных решений. В этом мы убедились при проведении соответствующих работ по исследованию алгоритмов обработки измерительной информации, получаемых в эталонах времени и частоты в 80-ые годы [9].
2. На рис. 4 приведены результаты адаптации коэффициентов фц, 012,013 при нулевых начальных значениях.
ÊL
<Al
к;
ISO S
Рис. 4. Процесс адаптации коэффициентов авторегрессии при нулевых начальных значениях
Процесс подстройки параметров носит примерно такой же характер, как и рассмотренный выше, только установившиеся (или почти установившиеся к конечному такту адаптации) значения коэффициентов 4h.ii 012, 01 э дальше, в смысле нормы разности, отстоят от оценок, полученных по «истинным» (т. е. сгенерированным) рядам.
Из выше сказанного следует, что процедура адаптации ни в коем случае не может подменять методику идентификации недоопределенных систем, детально рассмотренных в работе [3], а может служить только как дополнение (весьма важное) к алгоритму динамической обработки данных.
3. При работе с реальными объектами возможна «разладка» моделей. Маловероятно, что разладка может произойти для всех моделей одно-
временно. Поэтому при создании программного комплекса, реализующего рассмотренную выше процедуру адаптации моделей, необходимо предусмотреть возможность установки в начальное значение отдельных составляющих вектора шага р
(здесь мы рассматривали р как скалярную величину). Вопросы обнаружения разладки прогнозирующих моделей в данной статье не рассматривались, им посвящено достаточно большое число публикаций в специальных изданиях.
БИБЛИОГРАФИЧЕСКИЙ СПИСОК
1. Эльясберг П.Е. Определение движения по результатам измерений. М. : Наука, 1976. 416 с.
2. Ипполитов А.А, Хрусталев Ю.П. Субоптимальная фильтрация в системах с неполной матрицей наблюдений // Информационные и математические технологии в науке и управлении : тр. XV Байкал. Всерос. конф. Ч. I. Ируктск : ИСЭМ СО РАН. 2010. С. 174 182.
3. Хрусталев Ю.П., Акулов В.М., Ипполитов А.А., Курышева Л.Н. Обработка данных, полученных по результатам взаимных измерений вторичного эталона времени и частоты // Вестн. Иркут. гос. техн. ун-та. 2012. № 7. С. 22-29.
4. Ермольев Ю.М. Методы стохастического программирования. М. : Наука, 1976. 239 с.
5. Серышева И.А., Хрусталев Ю.П. Метод стохастического квазиградиента в задаче адаптации прогнозирующих моделей. // Вестн. Иркут. гос. техн. ун-та. 2013. № 12. С. 25-30.
6. Khrustalev Yu. P. Statistical and Dynamic Processing of Data Obtained when Handling Time and Frequency Standards // Measurement Techniques. 2004. V. 47. №. 6. P. 555-561.
7. Уидроу Б. Стационарные и нестационарные характеристики обучения адаптивных фильтров, использующих критерии минимума СКО / Б.Уидроу, Дж. М. Маккуи, М.Г. Ларимор, С.Р. Джонсон. // ТИИЭР. 1976. Т. 64. № 8. С. 37-51.
8. Демидович Б.П., Марон И.А. Основы вычислительной математики. СПб. : Лань, 2009. 664 с.
9. Хрусталев Ю.П., Спиридонова Е.В. Алгоритмы обработки измерительной информации, получаемой в процессе хранения единиц времени и частоты // Техника средств связи. Сер. Радиотехнические измерения. 1986. вып. 1. С. 58-72.
10.Химмельблау Прикладное нелинейное программирование. М. : Мир, 1975. 534 с.
11. Бокс Дж., Дженкинс Г. Анализ временных рядов, прогноз и управление / пер. с англ. В.Ф. Писаренко. М. : Мир, 1974.
