Построение динамических стохастических моделей, используемых при решении задач оценивания состояния групповых эталонов Текст научной статьи по специальности «Математика»

УДК 519.246.85+ 681.5.015.42 Хрусталёв Юрий Петрович,

к. т. н., доцент кафедры вычислительной техники ИрГТУ, тел.: 40-51-07, e-mail: [email protected]

Ипполитов Александр Александрович,

аспирант кафедры вычислительной техники ИрГТУ, тел.: 40-51-07, e-mail: [email protected]

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

Y.P Khrustalev, A.A. Ippolitow

BUILDING STOCHASTIC DYNAMICAL MODELS USED IN STATE ESTIMATION OF COLLECTIVE STANDARDS

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

Ключевые слова: недоопределённые системы, состояние, оценивание, статистика, модели, АРПСС.

Abstract. The problem of collective standards state estimation by measuring data obtained in process of their operation is considered. Collective standards, such as time and frequency standards, may be represented as linear dynamical systems with incomplete observation matrix. Findings obtained allow to estimate the state vector of collective standards in statical and dynamical data process modes.

Keywords: underdetermined system, state, estimation, statistics, models, ARIMA.

Задача оценивания состояния системы заключается в нахождении вектора её состояния Y по результатам выполненных в некоторый момент времени t наблюдений. Вектор наблюдений Z является функцией от вектора состояния, т. е.

Z(t ) = f (Y (t )). (1)

В линейных системах оператор L представляет собой матрицу H. В стационарных системах матрица H не зависит от времени. Задача оценивания состояния линейной системы может быть сформулирована следующим образом: по заданному вектору наблюдений Z = H ■ Y требуется найти оптимальную в некотором смысле оценку

Y вектора состояния Y. В качестве критерия оптимальности оценки обычно полагают минимум

суммы квадратов отклонений оценки вектора У от его истинного значения У для всех моментов времени tl, /=1,2.....N . В общем случае размерности матрицы H и вектора У не совпадают. Если размерность матрицы Н меньше размерности вектора У система является недоопределённой, т. е. системой с неполной матрицей наблюдений. В технических системах подобная ситуация имеет место довольно часто. Например, в групповых эталонах единиц физических величин один из элементов выбирается в качестве опорного. В определённые моменты времени производятся измерения разности физической величины, воспроизводимой каждым из элементов группового эталона с соответствующим значением опорного элемента. Схема измерений в эталоне приведена на рис. 1.

Рис. 1. Структурная схема измерений, выполняемых в групповом эталона

Нетрудно заметить, что в этом случае число измерений, выполненных в момент времени ti, на единицу меньше числа элементов группового эталона. Применение схемы измерений «каждый с каждым» не увеличивает размерность матрицы Н, т. к. получаемые при этом уравнения линейно зависимы от уравнений «основной схемы», т. е. схе-

У оп

n

Z gi (y*+ ^)

gi =

а„

(2)

а„

мы с опорным элементом. Примером подобных эталонов являются эталоны единиц времени и частоты. В эталонах времени измерения выполняются через определённые промежутки времени, их (эталоны) можно рассматривать как динамические системы. В качестве вектора состояний У (/) удобно выбирать отклонения частоты каждого из кван-тово-механических генераторов (водородных стандартов) от приписываемых им значений.

Задача нахождения оптимальных оценок вектора состояния по результатам наблюдения представляет собой задачу динамической фильтрации (если речь идёт об обработке данных в темпе их поступления). Решение этой задачи даётся системой рекуррентных уравнений, известных как фильтр Калмана [1].

Специфика групповых эталонов такова, что поведение каждого из его элементов не зависит от поведения других. Например, в эталонах времени и частоты каждый из водородных стандартов «развязан» от других по питанию, изолирован от механических нагрузок и т. п. Это позволяет достаточно просто произвести декомпозицию фильтра [2].

В полученном субоптимальном фильтре оценка состояния опорного элемента находится как линейная комбинация её составляющих: составляющая, найденная по результатам измерений, выполненных на данном такте измерений; составляющая, полученная на основе прогноза вектора состояния:

Х(t) = z(A)• fmp (t) + z(Б)•?(t) +

+ x{ В Ж t ) + #( t),

1, если факторы данного типа c

участвуют в формировании X (t),

(3)

где ^(c) =

0 - в противном случае,

где gi - вес каждого из элементов группового эталона; а2 - дисперсия прогноза i-й составляющей вектора состояния; уг - прогноз i-й составляющей, выполненный на предыдущем такте.

Т. к. результаты наблюдений представлены в виде временных рядов, то удобно в качестве прогнозирующих моделей выбрать именно модели временных рядов. Произведённая декомпозиция фильтра позволяет перейти от операций с многомерными временными рядами к работе с одномерными, что существенно упрощает процедуру обработки данных (на связь фильтров Калмана с моделями временных рядом указывал ещё академик В.С. Пугачёв [3]).

В общем случае временные ряды могут быть представлены в виде [4]

с = А, Б или В; / (1) - функция тренда (долговременная составляющая временного ряда - полиномиальный тренд); ф(1) - динамическая составляющая (динамический тренд); у(1:) - сезонная составляющая; ^) - случайная составляющая.

В выражении (3) ^) - случайная составляющая, значение которой в моменты времени ^, г = 1, 2,..., N, как правило, коррелированы. Это обстоятельство может быть использовано для построения стохастических динамических моделей -моделей авторегрессии - скользящего среднего (АРСС).

Эти модели наиболее широко используются при анализе временных рядов [5]. В частности, модели АРПСС успешно использовались при оценивании относительных отклонений частоты водородных генераторов [6]. Однако при разработке алгоритмов, опирающихся на прогнозы значений частоты, использовались модели АРСС, построенные по имеющимся временным рядам. Эти временные ряды - ряды относительных отклонений частот генераторов, входящих в эталон - были получены по результатам метрологического анализа выполненных измерений. При этом учитывалась вся доступная к моменту анализа информация, включающая в том числе и неформализованные процедуры. Методика построения моделей временных рядов, когда эти ряды имеются в наличии, широко известна [5]. Существуют реализации данной методики (методики Бокса - Дженкинса) в составе различных пакетов прикладных программ.

В нашем случае рассматривается более сложная задача: задача построения моделей АРСС по результатам косвенных измерений, когда в нашем распоряжении нет исходных временных рядов. Результаты измерений ^ представляют собой разности между относительными отклонениями частоты опорного и /-го элементов группового эталона в момент /. В этом случае непосредственное использование методики Бокса - Дженкинса невозможно. Поэтому предлагается использовать подход, основанный на минимизации функционала, представляющего собой сумму квадратов от-

1

i =1

клонений результатов измерений zi ) от их про-

гнозов:

N п—1

I=ЕЕк,—(?»(')—у^ (' ))]2,

(4)

]=1 '=1

где N - длина временных рядов, п - количество элементов группового эталона.

Т. к. прогнозы у. ) вычисляются на основе

моделей АРСС, то функционал I зависит от векторов параметров моделей в, т. е. I = I(в), где

РГ =[ф15 фР, ], Р - порядок авторе-

грессии (АР), q - порядок скользящего среднего (СС).

Параметры моделей должны удовлетворять условиям стационарности и обратимости, т. е. корни характеристических уравнений ф(в) = 0 и б(в) = 0 должны лежать вне единичного круга (В - оператор сдвига назад, ф(в) - оператор АР и 9(в) - оператор СС). Поскольку целевая функция - квадратичная, то при выборе начального приближения из области допустимых значений (например, % = 0,/ = 1, 2,..., р, в1 = 0,/ = 1, 2,..., </) оценки параметров ф, 91 могут быть найдены решением оптимизационной задачи без учёта указанных выше ограничений.

Проиллюстрируем использование такого подхода. При проведении машинного эксперимента для системы из трёх элементов (п = 3) для каждого элемента генерировались временные ряды, соответствующие моделям АРПСС размерности (1, 0, 0) с коэффициентами фф : фф = 0,3; фф = 0,5; фф = 0,7 - при значениях у\ = уф = уф = 0. Дисперсия шума, возбуждающе-

го систему, принималась равной сф = 0,1. Из этих рядов получались разностные ряды - ряды «измерений», которые затем обрабатывались с помощью вышеописанного алгоритма.

После достижения оптимального значения функционала I были получены оценки параметров авторегрессии:

фф = 0,36; фф = 0,46; фф = 0,67 В случае увеличения длины рядов наблюдений точность идентификации возрастала. Так, при увеличении выборки со 100 до 1000 наблюдений г-норма разности векторов истинных параметров моделей и их оценок уменьшилась ориентировочно в 4 раза [7]. Естественно, что с повышением сложности исходных моделей точность параметрической идентификации может снижаться.

Таким образом, при известной структуре моделей, т. е. при точно заданных порядках АР и СС (р и q) процесс подгонки параметров моделей минимизацией функционала I сходится. Однако точные сведения о порядках процессов АР и СС доступны только для сгенерированных данных и при работе с реальными рядами отсутствуют. Возникающая при этом проблема определения структуры моделей (Бокс и Дженкинс называют эту задачу идентификацией модели [5]) решается путём сравнения выборочных автокорелляцион-ной (АКФ) и частной автокорелляционной (ЧАКФ) функций с их теоретическими значениями. Несмотря на определённую субъективность понятий «похожести» и «неполной похожести» выборочных АКФ и ЧАКФ на теоретические функции моделей, у которых заданы порядки Р и q, метод даёт хорошие практические результаты и при определённых навыках исследователя позволяет достаточно точно определять структуру моделей АРСС. При этом авторы методики рекомендуют производить подгонку параметров для нескольких альтернативных моделей и выбирать наилучшую.

Поскольку в нашем распоряжении нет исходных временных рядов, а имеются только разностные временные ряды, такой подход к идентификации структуры моделей АРСС невозможен.

Можно предложить различные подходы к решению этой проблемы.

1. Простой перебор всех возможных вариантов структуры моделей с последующим сравнением. Исходя из того, что для реальных процессов порядки р и q весьма невелики (например, авторы [5] считают, что они не превосходят трёх, это предположение подтверждается на практике [6]), такой подход теоретически возможен. Однако даже при незначительных порядках АР и СС и при весьма ограниченном числе элементов, входящих в состав группового эталона, метод простого перебора требует очень больших затрат вычислительных ресурсов, поскольку необходимо многократно решать оптимизационную задачу. Компромиссный же вариант в виде использования одних и тех же значений порядка для всех моделей может не обеспечить соответствия выбранных порядков моделей «истинным», приводя к избыточности либо к снижению точности. Кроме того, при таком подходе возникает проблема разработки достаточно универсального критерия выбора «наилучшей» структуры модели. В результате проведённого моделирования такой подход признан малопродуктивным.

2. Результаты измерений представляют собой алгебраическую сумму двух процессов. Про-

Н =

1 0 0

-1

(5)

При числе элементов в групповом эталоне п ранг матрицы Н равен (п - 1). Решение уравнения

2 = Н • У (6)

в этом случае может быть найдено с помощью псевдообратной матрицы которая имеет вид

[7]

Н+= Нт •(Н• Нт)-1. (7)

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

3. В случае когда имеются априорные сведения об исследуемых процессах, возможно использовать их для структурной идентификации моделей. В частности, в 80-е годы был проведён большой объём работ по построению моделей АРСС для рядов относительных отклонений частот водородных генераторов, эксплуатируемых в составе эталона-копии Государственного эталона времени и частоты СССР [6]. При этом следует учитывать, что состав аппаратуры эталона с тех пор существенно изменился (менялись как стандарты частоты, так и измерительный комплекс). Кроме того, даже для того состава эталона-копии модели временных рядов для различных водородных генераторов (стандартов частоты) существенно отличаются друг от друга. Поэтому и данный подход вряд ли может быть признан продуктивным.

4. Можно построить по результатам измерения некоторые приближённые («грубые») оценки рядов =1, 2, ..., N . Затем по полученным временным рядам, используя методику Бокса-Дженкинса, построить модели АРСС. Оценки параметров этих моделей, в том числе и оценки порядков р и q для каждой из моделей, а при необходимости и оценки коэффициентов при членах АР и СС принимаются затем в качестве начального приближения при решении оптимизационной задачи.

Приближённые оценки временных рядов у (^ ) можно найти методом наименьших квадратов. Матрица наблюдений Н для каждого из моментов времени имеет для групповых эталонов следующий вид (без потери общности результатов будем считать опорным 1-й элемент эталона): 1 -1 0 ••• 0 1 0 -1 ••• 0

Таким образом, получаем:

Н+ =-• п

1

-(п -1)

1

1

-(1,-1)

1

(8)

Тот же самый результат может быть получен с помощью метода Гревиля [8].

Решение уравнения (6) даётся выражением

У = Н+ 2 (9)

для каждого из моментов (г = \, 2, ..., Л^)

В частности, для опорного элемента имеем

1 П-1

У = - •£ ^. (Ю)

п ,=1

Оценки уг-, (г = 2, 3,..., п) могут быть найдены из соотношения

= У - Уг (11)

подстановкой вместо величины у1 её оценки (10) либо непосредственно из формулы (9) (как упоминалось выше, шумами измерений в данном случае пренебрегаем, что вполне допустимо при обработке данных на суточных интервалах). Погрешность оценивания значения опорного элемента при этом равна:

1 п

А = 1 •£ У,. (12)

П 1=1

Или, если использовать взвешенное среднее:

п

А 2 =Е ЯУ. (13)

,=1

Естественно, что условия равенства выражений (12) и (13) нулю выполняются далеко не всегда, поэтому при применении такого подхода временные ряды, по которым строятся модели АРСС, содержат погрешности. Анализ этих погрешностей не входит в рамки данной работы. В целом можно считать вполне приемлемым применение моделей АРСС, полученных на основе изложенного подхода, в качестве начального приближения при решении оптимизационной задачи:

«Найти значение расширенного вектора Р = [р!, Р2,..., Рп], доставляющее минимум функционалу (4)». р1 =[ф1,ф2,...,фр;в;,е2,.,е;] -

вектор параметров модели АРСС /-го элемента группового эталона.

Чтобы уменьшить влияние погрешностей на качество построенных по предлагаемой методике моделей, можно решать оптимизационную задачу при разных начальных приближениях. Это соответствует рекомендациям Бокса и Дженкинса по подгонке параметров нескольких альтернативных моделей.

Продемонстрируем использование данной методики структурной идентификации моделей АРСС с помощью машинного эксперимента. Сгенерируем ряды «наблюдений» по некоторым заданным моделям АРСС (наложив на них псевдослучайный шум). Затем методом наименьших квадратов построим ряды «грубых» оценок. Анализируя их АКФ и ЧАКФ (пример АКФ и ЧАКФ приведён на рис. 2-3), определим порядки процессов АР и СС. Выполним процедуру подгонки моделей и построим оценки исходных рядов. Количество исходных рядов в эксперименте - 4, количество полученных из них разностных рядов - 3, длина рядов - 300. Параметры моделей и полученные оценки параметров процессов АР и СС приведены в табл. 1.

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

В результате анализа коррелограмм, для ряда 1 можно предположить структуру с р = 0; q = 1, т. е. (0, 0, 1). Для ряда 2 следует предпочесть мо-

дель вида (0, 0, 1). Как альтернативную можно рассмотреть и (1, 0, 0). Меньшее число членов, чем в исходной модели, здесь может быть объяснено относительной малостью соответствующего коэффициента исходного ряда и, соответственно, малостью его вклада в процесс СС. Для ряда 3 следует выбрать структуру модели (3, 0, 0), а для ряда 4 - (2, 0, 0).

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

Используя эти результаты (табл. 1), можно перейти к процедуре оценивания. Как можно увидеть из таблицы, оценки коэффициентов авторегрессии и скользящего среднего при заданных порядках процессов достаточно близки к истинным.

Приведём график (40 значений) рядов «грубых» оценок и рядов уточнённых оценок для опорного элемента, построенных по предложенной методике (рис. 4). Оценки для остальных элементов вычисляются на основе оценки опорного элемента и поэтому не требуют отдельного рассмотрения.

Рис. 2. АКФ, ЧАКФ исходного ряда №3

Рис. 3. АКФ, ЧАКФ рядов «грубых» оценок для ряда №»3

Т а б л и ц а 1

Номер ряда 1 2 3 4

Порядок авторегрессии Исходный Р1 = 0 Р2 = 0 Рз = 3 Р4 = 2

Оценка Р1 = 0 Р2 = 0 Р3 = 3 Р4 = 2

Коэффициенты авторегрессии Исходные - - Ф1 = 0,5 ф2 = -0,3 ф3 = -0,2 Ф1 = 0,5 ф2 = -0,4

Оценки - - ф! = 0,275 ф2 = -0,255 ф3 = -0,217 ф1=0,531 ф2= -0,426

Порядок скользящего среднего Исходный 41 = 1 42 = 2 43 = 0 44 = 0

Оценка 41 = 1 42 = 1 43 = 0 44 = 0

Коэффициенты скользящего среднего Исходные 01 = 0,35 01=0,2 02=-0,1 - -

Оценки 0! = 0,237 01=0,159 02=0 - -

с. к. о. псевдослучайного шума Ох = 0,01 а2 = 0,02 а3 = 0,01 а4 = 0,02

Рис. 4. Сопоставление исходного ряда (ромбы), ряда оценок среднего (круги) и ряда уточнённых оценок (квадраты) для опорного элемента

Из графика видно, что уточнённые оценки лучше описывают исходный временной ряд, чем оценки, найденные простым применением метода наименьших квадратов, - «грубые оценки». Средний квадрат отклонений уточнённых оценок (найденный по трёмстам точкам) от истинных значений ряда равен 0,0029. Соответствующая величина для «грубых» оценок равна 0,0044, т. е. в полтора раза больше.

Мы рассмотрели методику построения моделей АРСС для случая, когда в формировании временного ряда участвует только составляющая

^(ь). Чтобы воспользоваться предложенной методикой при наличии трендов, необходимо найти оценки / (1), ф(1) и ) (см. выражение 3)

и учесть их вклад в ряды ().

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

БИБЛИОГРАФИЧЕСКИЙ СПИСОК

1. Ли. Р Оптимальные оценки, определение характеристик и управление. М. : Наука, 1966. 176 с.

2. Ипполитов А. А., Хрусталёв Ю. П. Субоптимальная фильтрация в системах с неполной матрицей наблюдений. // Информационные и математические технологии в науке и управле-

нии : тр. XV Байкальской Всерос. конф. Ч. 1. Иркутск : ИСЭМ СО РАН, 2010. С. 174-182.

3. Пугачев В. С. Теория вероятностей и математическая статистика : учеб. пособие. 2-е изд., испр. и доп. М. : Физматлит, 2002. 496 с.

4. Айвазян С. А. Прикладная статистика. Основы эконометрики. Т. 1 : Теория вероятностей и прикладная статистика. / С. А. Айвазян, В. С. Мхитарян. М. : Юнити-Дана, 2001. 656 с.

5. Бокс Д., Дженкинс Г. Анализ временных рядов. Прогноз и управление. Вып. 1. М. : Мир, 1974. 406 с.

6. Хрусталёв Ю. П., Спиридонова Е. В. Алгоритмы обработки измерительной информации, по-

лучаемой в процессе хранения единиц времени и частоты // Техника средств связи. Сер.: Радиотехнические измерения. 1986. Вып. 1. С.58-72.

7. Хрусталёв Ю. П., Овечкина А. А., Щербаков Е. В. Построение моделей многомерных временных рядов по результатам наблюдений в динамических системах. // Методы исследований и моделирования технических, социальных и природных систем : сб. науч. тр. Новосибирск : Наука, 2003. С. 293-307.

8. Гантмахер Ф. Р. Теория матриц. М. : Наука, 1966. 576 с.

УДК 519.685 Богданова Вера Геннадьевна,

к. т. н., с. н. с., Федеральное государственное бюджетное учреждение науки Институт динамики систем и теории управления Сибирского отделения Российской академии наук (ИДСТУ СО РАН), лаборатория параллельных и распределенных вычислительных систем, e-mail: [email protected]

Горский Сергей Алексеевич,

к. т. н., н. с., Федеральное государственное бюджетное учреждение науки Институт динамики систем и теории управления Сибирского отделения Российской академии наук (ИДСТУ СО РАН) , лаборатория параллельных и распределенных вычислительных систем, тел. 89149230010, e-mail: [email protected]

ТЕХНОЛОГИЯ ПАРАЛЛЕЛЬНОГО РЕШЕНИЯ НЕЛИНЕЙНЫХ СИСТЕМ БУЛЕВЫХ УРАВНЕНИЙ НА ВЫЧИСЛИТЕЛЬНОМ КЛАСТЕРЕ

V. G. Bogdanova, S.A. Gorskiy

THE TECHNOLOGY OF PARALLEL SOLUTION OF NONLINEAR SYSTEMS OF BOOLEAN EQUATIONS ON A COMPUTING CLUSTER

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

Ключевые слова: параллельные вычисления, булево моделирование, решение булевых уравнений.

Abstract. Enterprise principles, base functions and structural components of toolkit for automation of parallel computing of Boolean equations on a homogeneous computing SMP cluster is considered. In particular, the parallel solution of Boolean equations based on the hybrid approach is discussed.

Keywords: parallel computing, Boolean modeling, Boolean equations solving.

Введение

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