УДК 519.246.8
МЕТОДИКА СТРУКТУРНО-ПАРАМЕТРИЧЕСКОЙ ИДЕНТИФИКАЦИИ СИСТЕМЫ ВРЕМЕННЫХ РЯДОВ
© 2013 Ю.Е. Кувайскова
Ульяновский государственный технический университет
Поступила в редакцию 26.09.2013
Описывается методика структурно-параметрической идентификации системы взаимосвязанных временных рядов. Методика основана на подходе адаптивного динамического регрессионного моделирования, предусматривающего при построении математических моделей временных рядов последовательную адаптацию модели к возможным нарушениям основных предположений регрессионного анализа. При этом построенная модель системы временных рядов адекватна реальной ситуации и позволяет прогнозировать с большей точностью будущее состояние процесса.
Ключевые слова: система временных рядов, адаптивное динамическое регрессионное моделирование, прогнозирование.
Одной из важных задач во многих отраслях промышленности является задача прогнозирования будущего состояния технического объекта. В большинстве случаев состояние объектов в различных технических приложениях характеризуется набором значений параметров, фиксируемых через определенные моменты времени и образующих систему временных рядов. При этом часто наблюдается высокая степень корреляционной зависимости между временными рядами значений параметров технического объекта. При моделировании динамики временных рядов нередко сталкиваются с такими нарушениями схем регрессионного анализа, как высокая степень автокорреляционной зависимости между последующими и предшествующими членами временного ряда, ненормальность распределения остатков и другие. В связи с этим возникает задача комплексной обработки системы взаимосвязанных временных рядов с обязательной проверкой соблюдения предположений нормальной схемы Гаусса-Маркова и последующей адаптацией к их нарушениям.
Для решения данной задачи разработана методика структурно-параметрической идентификации системы взаимосвязанных временных рядов на основе методологии адаптивного динамического регрессионного моделирования [1, 2]. Подход адаптивного динамического регрессионного моделирования представляет собой реализацию метода многоэтапной структурно-параметрической идентификации. На каждом этапе проводятся построение и анализ соответствующей компоненты модели временного ряда, оценка точности аппроксимации и прогнозирования,
Кувайскова Юлия Евгеньевна, кандидат технических наук, доцент кафедры «Прикладная математика и информатика». E-mail: [email protected]
диагностика свойств остатков и при необходимости - адаптация к нарушениям предположений о характеристиках этих свойств.
Рассмотрим технический объект, состояние которого характеризуется параметрами, значения которых регистрируются через определенные промежутки времени и образуют систему временных рядов, которую можно записать в виде У(г) = МО, ЖХ ■■■, Ж)} , где параметр г указывает на момент времени, в который зафиксировано значение, или номер наблюдения.
Пусть г = 1, 2, ■..,п (т.е. количество наблюдений равно п), рассмотрим N временных рядов
у1 (г), у2 (г), ■.., у* (г).
В соответствии с подходом адаптивного динамического регрессионного моделирования одномерный временной ряд у' (г) , наблюдаемый в равноотстоящие моменты времени г1, г2, ■.., гп представляется в виде:
У'(г) = Г' (г) + ё (г) + У (г) + (г), (1)
при этом функция Г'(г) - неслучайная (долговременная) функция тренда соответствующего ряда; ё'(г) - неслучайная периодическая функция - совместная гармоническая составляющая ряда; у' (г) - случайная с элементами регулярности функция - векторная авторегрессия; е'(г) - нерегулярная компонента (случайная величина, ошибка).
На первом этапе для системы временных рядов выполняется разведочный анализ (фрактальный, мультифрактальный) для выявления степени регулярности каждого ряда системы [3].
При выявлении заметной трендоустойчиво-сти для соответствующих рядов выделяется функция тренда Г(г), которою обычно приближают полиномом достаточно низкой степени, также могут быть использованы и различные нелиней-
ные по оцениваемым параметрам выражения. Оптимальная функция трендовой составляющей ищется по критерию минимума внешнего сред-неквадратического отклонения (СКО) <А .
Точность аппроксимации < (внутреннее СКО) рассчитывается по формуле:
а
1
Ё(y, - y,)2/(n - p),
(2)
i=1
°А=, Ё (Д,-A)2/(k -1)'
(3)
где А, = y{ — yi, i = 1, k (k - количество элементов контрольной выборки),
А = -
Ё А,
=1_
k
(4)
у которых на данном периоде Тк имеются значимые гармоники.
Совместную гармоническую составляющую можно представить в виде:
где п - количество наблюдений,р - число слагаемых в модели, у1 - результат ¿-го наблюдения, у1 - прогнозируемое значение по построенной модели.
Для вычисления внешних мер (внешнего СКО <А и др.) [1] исходная выборка данных делится на две части - обучающую и контрольную. По обучающей выборке строится модель, по контрольной - внешние меры, характеризующие прогностические свойства модели.
Точность прогнозирования <А вычисляется по формуле:
yi - наблюдаемое значение отклика на контрольном интервале, yi - его прогноз - значения, вычисляемые по комплексной модели.
При этом для оценивания применяется линейный или нелинейный метод наименьших квадратов. Полученная модель для одного ряда может быть представлена в виде:
У' (t) = fi (t) + £_ trend' (t), (5)
где y' (t) - фактические значения ¿-го ряда; f' (t) - функция трендовой составляющей для ¿-го ряда, е' trend(t) - случайная величина, характеризующая отклонения реального значения ряда от теоретического, найденного по уравнению регрессии.
Далее из каждого ряда системы выделяются остатки е _ trend1 (t):
е _ trend'(t) = y' (t) — f'(t). (6)
Дальнейшее сглаживание временных рядов производится методами гармонического анализа, основанного на использовании формул Фурье.
Для каждого временного ряда методом пошаговой регрессии находятся значимые гармоники. Затем последовательно ищутся совместные амплитуды и фазы только для тех временных рядов,
n/2
(
g (t) = Ё A ■ sin
k=1
2nt
Л
T
V 1k
(7)
Ё A,
где t = 1, 2,n, a := J=L
j ,k
T
Коэффициенты
Л]к, р]к и = 1, 2,..., М, к = 1, 2, ..., п)
находим для каждого из N временных рядов по формулам:
Ajk = V w'k+v %
Vjk = arctg (j)
j wjk
w
jk
= Ё е _ trend
( 2nd ^
jd
■ Sin
T,
k у
jk
= Ёе_trendjd ■ cos
( 2nd Л
d =1
T
V 1k у
(8)
(9)
(10)
t = n Tk k.
Далее находятся остатки для каждого временного ряда:
е_g' (t) = е_trend (t) — g' (t). (11)
После выделения регулярных составляющих целью анализа системы временных рядов является моделирование остатков е g' (t) случайной с элементами регулярности функцией y(t) . Функция y/(t) для системы взаимосвязанных временных рядов может быть представлена в виде модели векторной авторегрессии.
Модель векторной авторегрессии - это система уравнений, каждое из которых представляет собой модель авторегрессии и распределенного лага. Пусть x'(t) ' = 1..N - ¿-й временной ряд. Модель векторной авторегресии порядкаp - модель для ¿-го временного ряда будет иметь вид
x'' (t) = a0 + Ё j (t — 1) +
j=1
Ё a2jXJ (t — 2) + ... +
j=1
N
+ Ё a'pjXJ (t — p ) + е (t)
j=1
'=1
d =1
Матричная запись модели векторной авторегрессии порядка p для системы временных рядов Х(У) = {(У), x2(У), ..., хм(У)} запишется в виде
Х(У) = а0 + Л1Х(У -1) + Л2Х(У - 2) +... +
р
+ АрХ (У - р) + е(У) = а0 + X ЛтХ(У - т) + е(У)
т=1
(13)
(1 х'(у - 1) 1 х'(У - 1)
Х Ч У - р) хЧУ - р)
ХИ (У - 1) ... ХИ (У - р )
Xя (у - 1)
1 X'(У - 1) ... X'(У - р ) ... Xя (У - 1) ...
Xя (У - р) Xя (у - р)у
,(15)
единичный столбец соответствует нулевому коэффициенту.
Y- вектор последних значений всех рядов:
(X1 (У)...Xя(У))Т . (16)
В результате операций над матрицами, получаем матрицу коэффициентов A, имеющую (1+^р) строк и N столбцов.
Г а0
V ая - р
а,
я
а
(17)
я ■ р у
Каждый столбец здесь соответствует коэффициентам для ¿-ого уравнения модели.
Остатки £г g(у) после вычитания из исходного ряда регулярной составляющей представляются в виде модели векторной авторегрессии:
я
£_ gi (У) = а0 + X а-^у' (У - 1) +
' =1
(18)
В итоге получаем для каждого из N временных рядов комплексную модель следующего вида:
у' (У) = /' (У) + gi (У) + а 0 + + ¿а '^у' (У - 1) + ¿а 2'У' (У - 2) + ... +
(20)
'=1
'=1
р
где Лт - матрицы элементов ат' ' .
Коэффициенты уравнений модели можно получить методом наименьших квадратов, пользуясь следующей формулой:
Л = (ХТХ)-1 ХТУ. (14)
Здесь X - это матрица значений рядов, представленная в виде:
+ X а2'У' (У - 2) + ... +
'=1
я
+ X a'pjУJ (у - р) + '=1
+ ¿,а'ье_ gi (У - к) +е' (У)
Оптимальная модель векторной авторегрессии (18) ищется по критерию минимума внешнего СКО <УА , рассчитанного по формуле (3). Далее из каждого ряда системы выделяются остатки:
( я я я р Л
* (У) = ^ (У) - и + Xа-'У' (У-1) + ¿а2 у (У - 2) +...+ ¿а' (У - р) + ¿агье_ ^ (У - к) I.
V '=1 '=1 '=1 к=1 У
+ X а'р'У' (У - р) + X а'и£ _ g' (У - к) + е' (У)
'=1 к=1
Для построенных моделей по формулам (2) и (3) подсчитывается точность аппроксимации < и прогнозирования <А.
После структурно-параметрической идентификации системы временных рядов проверяется соблюдение условий применения "регрессионного анализа-метода наименьших квадратов". Если основные предположения соблюдаются, построенная комплексная модель системы временных рядов может быть использована для прогнозирования.
При несоблюдении предположений проводится соответствующая процедура адаптации к нарушению данного предположения [1].
Ниже представлены результаты анализа, моделирования и прогнозирования временных рядов вибраций гидроагрегата в режиме работы в сети. Наблюдения в режиме работы агрегата в сети длились 1200 секунд, показания снимались 2 раза в секунду, всего получено 2400 наблюдений по каждому из восьми каналов. После усреднения каждые 10 секунд получилось 120 наблюдений.
Применяя методику структурно-параметрической идентификации системы временных рядов для моделирования измерений вибраций по восьми каналам, получены следующие модели:
у1 (У) =641946+ 0,025 ■ У +1,022 ■ 5ш(- 26,68 ■ У + 65,093) + + 0,001 у2(У-1 >0,001 у3(У4) - 0,002-у4(У-1) + + 0,001- у 5(У-1) + 0,0004- у 6(У-1 >0,002- у 7(У-1) + в1 (У),
у2(У) =698,017 - 0,003 ■ У +1,022 ■ 5т(12,493 ■ У - 22,256) -
- 0,0002- у1(У-1 >0,006- у3(У-1) + 0,007- у4(У-1) -
- 0,004 - у5(У-1 >0,001- у6(У-1) + е2(У),
уЪ(У)=225,93 + 0,997 - 51п(0,0001 - У +11,456)- 0,317 - у1(У-1) + 0,002 - у4(У-1) + 0,0001 - у5(У-1) -
- 0,001 - у6(У-1) + 0,011 - у7(У-1) + 0,005 - у8(У-1) + еъ(г),
у4(У)=44,752- 0,109- У +1,021- 81п(236654- У +170546) + + 0,001- у1(У-1) + 0,0002- у2(У-1) + 0,008- у3(У-1) -
- 0,001- у5(У-1) + 0,002- у6(У-1 )-0,014- у7(У-1) -
- 0,007- у8(У-1) + е4(У),
к =1
y5(t) =29,174 + 0,008 • t + 1,023 • sin(-1,299 • t + 0,041) -
- 0,001 • y1 (t-1 >0,0001 • y 2(t-1) + 0,0001 • y 3(t-1) -
- 0,0003 • y6(t-1) + 0,0002 • y8(t-1) + e5(t),
y6(t) = 3,943 • t + 1,014 • sin(- 0,191 • t + 2,518) + + 0,109 • y Yt-1)-0,001 • y2(t-1) - 0,199 • y 3(t-1) -
- 0,024 • y4(t-1)-0,002 • y5(t-1)-0,129 • y7(t-1) + e6(t),
y 1(t) = 26,049 + 1,001 • sin(0,040 • t - 0,011)- 0,001 • y '(t-1) + 0,001 • y 3(t-1) + 0,0002 • y 4(t-1) -
- 0,001 • y 5(t-1 )-0,001 • y 6(t-1) + 0,004 • y8(t-1) + e7(t),
y 8(t)=32,155 + 1,022 • sin(- 3197,496 • t + 734,684)- 0,0004 • y1 (t-1) + 0,016 • y3(t-1) - 0,008 • y4(t-1) + + 0,0003 • y 5(t-1) + 0,0004 • y 6(t-1 >0,004 • y 7(t-1) + e 8(t).
На рис. 1 а показаны графики исходных данных (сплошная линия) и аппроксимации по соответствующей модели (штриховая линия), а также графики (Рис. 1 б) предсказанных (штриховая линия) и исходных (сплошная линия) значений временных рядов для данных вибраций по первому каналу.
Для сравнительного анализа эффективности разработанной методики моделирования системы временных рядов использовалась методика построения моделей авторегрессии проинтегрированного скользящего среднего (АРПСС), реализованная в пакете Statistica [4].
Построенные модели АРПСС в пакете Statistica, представлены ниже.
у1 (/)=643,6562+0,2151у '(/-1), у 2(/)=697,9956-0,3888 • е 2(/-1 >0,1814 ■ е 2(/-2), у ъ(/) =2116052+0,51052 ■ у 3(/-1), у 4(/)=39,69158-0,10275 ■ е 4(/-1)+0,20859 ■ е4 (/- 2), у5(/)=28,76911-0,20671 ■ е5(/-1 )+0,12513 ■ е5(/-2), у6(/)=-42,3635+0,8185 ■ у 6(/-1 )+0,1815 ■ у6(/-2) -- 0,0973 ■ е6(/-1)+0,5629 ■ е6(/-2), у7 (/)=25,28663+0,96942 ■ у7 (/-1)+ + 0,61375 ■ е7(/-1)+0,22012 ■ е7(/-2), у8(/)=31,82976-0,06093 ■ е8(/-1).
Результаты сравнения по точности прогнозирования (внешнее СКО) двух методик для системы временных рядов, полученных по восьми каналам измерения вибраций в режиме работы в сети гидроагрегата, приведены в табл. 1.
Таким образом, совместное описание характеристик технического объекта с применением методики многоэтапной структурно-параметрической идентификации с адаптацией к нарушениям основных предположений регрессионного анализа повышает точность прогнозирования по сравнению с точностью при использовании одномерных методов (методика АРПСС), что позволяет принимать с определенным запасом времени управленческие решения, направленные на предупреждение опасной ситуации.
Рис. 1. Моделирование (а) и прогнозирование временного ряда (б)
Таблица 1. Оценка точности прогнозирования для системы временных рядов при использовании двух методик
~"—--...__^Номер канала Мето дик а " Среднеквадратическое отклонение прогнозирования (сд)
1 2 3 4 5 6 7 8
Методика АРПСС
^аЬяНса) 4,99 3,3 2 2,68 15,91 8,21 878,55 6,01 2,72
Методика структурно-
параметрическои иде нт ифи к аци и систем ы
временных рядов 3,56 2,57 2,25 14,17 6,12 667,124 3,60 2,20
Исследование выполнено при поддержке Министерства образования и науки Российской Федерации, соглашение 14.B37.21.0672 (федеральная целевая программа "Научные и научно-педагогические кадры инновационной России").
СПИСОК ЛИТЕРАТУРЫ
1. Валеев С. Г. Регрессионное моделирование при обработке наблюдений. М.: Наука, 1991. 272 с.
2. Валеев С. Г., Кувайскова Ю. Е. Программное обеспечение обработки временных рядов техногенных характеристик // Обозрение прикладной и промышленной математики. 2009. Т. 16. Вып. 6. С. 1037-1038.
3. Валеев С. Г., Кувайскова Ю. Е., Губайдуллина С. А. Применение мультифрактального анализа при описании временных рядов в технике и экономике // Вестник Ульяновского государственного технического университета. 2008. №. 2. C. 23-27.
4. Боровиков В.П., Ивченко Г.И. Прогнозирование в системе STATISTICA в среде Windows. М.: Финансы и статистика, 1999. 384 с.
METHOD OF STRUCTURAL PARAMETRIC IDENTIFICATION OF TIME SERIES SYSTEM
© 2013 Ju.E. Kuvayskova
Ulyanovsk State Technical University
Method of structural parametric identification of relative time series system is described. The method based on procedure of adaptive dynamic regression modeling, providing at construction of mathematical models of time series consecutive adaptation of model to possible infringements of the basic assumptions of regression analysis. Thus constructed model of time series system is adequate to a real situation and allows predicting a process condition with greater exactness
Keywords: time series system, adaptive dynamic regression modeling, predicting.
Julia Kuvayskova, Associate Professor at the Applied Mathematics and Computer Science Department. E-mail: [email protected]