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

УДК 681.513.7:519.7

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

БОДЯНСКИЙ Е.В., ВОРОБЬЕВ С.А., КОСТЮК О.В., ЛЮБЧИК Л. М.

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

Задача обработки нестационарных стохастических последовательностей с тренд-сезонной компонентой достаточно часто встречается на практике и, прежде всего, в экономическом прогнозировании и технической диагностике [1-9]. При этом предполагается, что анализируемый ряд может быть представлен в виде модели

Р ~ m

Ук = Xdjk1 + Xcos а к +

,=0 j=1

+ ~j sin0jk) +j~k,

или через оператор сдвига назад z-1:

m

(1 - z_1)p П(1 - 2cos® jZ-1 +

j=1

+ z ~2)Ук =<~k ■

Здесь p и m — порядок полиномиальной компоненты и количество гармоник в последовательности у к соответственно; dt, a j , bj — коэффициенты модели, в общем случае неизвестные; О < а j = 2nfj T0 < п — неизвестные частоты гармонических компонент, при этом а, Фа j, i Ф j; T0 — период квантования; к = 1,2,. N — дискретное время; £к — стохастическая компонента типа

белого шума с нулевым математическим ожиданием и ограниченным вторым моментом. В общем случае модель должна иметь несмещённую структуру (значение m известно заранее), что, однако, не является обременительным ограничением.

Для оценивания неизвестных параметров модели (1) в [1-9] предлагается использовать метод наимень-

(1)

(2)

тих квадратов, при этом частоты а j полагаются

известными, поскольку входят в описание ряда нелинейно и их оценивание сопряжено со значительными трудностями. Однако в ряде случаев именно частоты представляют особый интерес, в связи с чем необходимо предусмотреть возможность их оценивания параллельно с коэффициентами dt, a j , bj .

В [10-12] предложена многоэтапная схема оценивания параметров полигармонических сигналов, которая для нашего случая может быть представлена в следующем виде:

1. Исключение полиномиальной компоненты путем взятия (р +1) -й разности последовательности

~ ~ р у ,

у к или вычитание из у к компоненты X dik , где

i=0

d. — оценки, получаемые с помощью метода наименьших квадратов.

2. Оценка параметров модели

m

П (1 - 2cosjjZ _1 + z ~2)Ук =£к

j=1

с помощью метода наименьших квадратов.

3. Восстановление оценки частот а j путем нахождения корней полинома m -й степени.

4. Нахождение оценок а j, b j с помощью метода

наименьших квадратов и пересчет их в оценки а j , b j .

5. Дигрирование центрированного сигнала у к и

получение отфильтрованной оценки ряда у к.

В качестве иллюстрации данного подхода рассмотрим несколько простых примеров.

Пусть в (1) р = 0, т. е. анализируемая последовательность у к колеблется относительно среднего

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

m

у к = Z (aj cosаJk +bjsin аJk) + £к> (3)

j=1

или

m

П (1 - 2cosj jz -1 + Д2)Ук =€к , (4)

j=1

іде у к = у к - У к-1, £к =<ук -<ук-1. При этом несложно видеть, что существует однозначная связь между коэффициентами

a j = a j (1 - cos а j) + bj sin а j,

* у у (5)

bj = bj (1 - cos®j) - aj sin coj.

74

РИ, 1998, № 3

Оценивание параметров модели (3), (4) с помощью метода наименьших квадратов осложняется тем,

что шум уже не является белым, что может

привести к смещённым (ридж-) оценкам [11], хотя в [12] и утверждается, что такое смещение не возникает. Так или иначе, в ряде случаев удобнее из yk вычесть среднее

- у 1 N ~

d0 = У = —I yk (6)

N k=1

и перейти к центрированному сигналу yk с белым

шумом £k. Этот процесс можно организовать в реальном времени, вычисляя на каждом такте среднее

~ у 1 ~ k —1 ~

do,k = yk =7 yk +-Т- yk-1 (7)

kk

и центрированное значение yk = у, — у, . Если

средний уровень d0 варьирует во времени, то вместо (7) можно использовать экспоненциальное среднее

d0,k = ~k = ayk + (1 — a)yk—1 , (8)

где 0 <a < 1 — параметр сглаживания [2, 7]. Далее переходим к пункту 2 многоэтапной схемы оценивания.

Пусть тестовый временной ряд генерируется с помощью следующей модели:

2

yk = d0 + I sin 2nfjTok + %k ■

J =1

Здесь do = 2,4 ; /1 = 50 ; /2 = 120 ; To = 1 —

квант машинного времени; — последовательность случайных величин с нулевым средним и

среднеквадратическим отклонением о^ = 2 , т. е.

генерируемый гармонический ряд сильно искажен действием помехи и смещен относительно центра на величину 2,4.

Ряд yk показан на рис. 1. На рис. 2 приведена

первая разность ряда yk. Очевидно, что первая разность колеблется относительно нуля. Таким образом, последовательность yk = yk — yk—1 является центрированной. На рис. 3 также представлены среднее ряда yk, вычисленное по формуле (7), и

центрированный ряд yk = yk — d0k . Из-за действия интенсивной помехи среднее ряда j^k варьирует во времени.

Пусть теперь в (3) m = 1, т. е. в сигнале присутствует лишь одна гармоника с частотой а> . Тогда (4) можно переписать в виде

(1 — 2 cos т — + z ~2) yk =^k, (9)

(1 — в2 z-1 + z-2) yk =4, (10)

о

Рис.1. Временной ряд yk

yk = yk — yk—1

с

Рис.3. Среднее d0k и центрированный ряд

yk = yk — d0,k

yk =в2 у,—1 — yk—2 +£k (11)

и ввести критерий идентификации

N

J1 = I ((yk + yk—2) ~ P2yk—i) , (12)

k=3

минимизация которого ведёт к оценке

N

I( у, + Уk—2 ) У k—1

в = —-----N---------■ (13)

21 yk—1

k=3

Заметим, что если у, есть p -я разность

процесса у,, то нижний индекс суммирования в

(12) начинается с k = p + 3 .

Далее несложно найти оценку частоты

т = arccos в, (14)

отфильтрованный сигнал

yk, N = 2вуk—1 — yk—2 (15)

РИ, 1998, № 3

75

и прогноз на один шаг

У N+1 = 2вуЫ - Уы-1- (16)

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

N

J2 = Z (Ук - a cos Ок - b sin Ок)2 =

к=3

N

f

= Z (Ук -(a b)

к=3

cos

сок Л

V sin Ок j

(17)

и получить оценки коэффициентов в виде

(а r

Vb J

Z cos2 Ок Z cos сок sin Ок

Z cos^эk sinc^ Z sin2 Ок

^Z cos^эkyk^

Z sin Окук

Y

(18)

2

)

x

x

Задача оценивания параметров о , а и b может быть решена в реальном времени, для чего целесообразно воспользоваться экспоненциально взвешенной процедурой стохастической аппроксимации [1315], обеспечивающей компромисс между фильтрующими и следящими свойствами процесса идентификации. В этом случае, если по N наблюдениям были

получены оценки О N , а N , bN , то с приходом

N +1 -го наблюдения производится уточнение согласно рекуррентной схеме

вN+1 = Pn + rN+1 (yN+1 + yN-1 -- Pn 2yN )2yN,

rN+1 = Yn + 4УІ> 0 — Y —1, (19)

О N+1 = arccos PN+b

а

N+1

V bN+1 J

Л ra..л

N

V bN J

+ RN+1 (yN+1 аы X

x cos О N+1(N +1) - bN sin О N+1(N +1)) x

cos

О N+1 (N + 1) ^

(20)

sin О N+1 (N +1)

rn+1 =7rn +1 0 — Y — 1-

Несложно видеть, что при у = 0 (20) превращается в градиентную процедуру с единичным шагом, а при у = 1 - в алгоритм Кифера-Вольфовица.

Заметим, что в многоэтапной схеме результаты каждого этапа влияют на точность последующего,

поэтому неточность в определении оценки в? влечёт за собой накопление погрешностей в оценках О , а,

I). Поэтому если для решения задачи фильтрации и прогнозирования достаточно лишь в (см. (15), (16)),

то этим и следует ограничиться.

Рассмотрим ещё один численный пример. В этот раз временной ряд генерируется моделью вида

Ук = 0,93sin 2 л/7ок + %к •

Здесь f = 50, £к — последовательность случайных величин с нулевым средним и дисперсией, равной 0,4.

Ряд у к и его прогноз у к приведены на рис. 4.

На рис. 5 показан средний квадрат ошибки прогнозирования.

Рис.4. Ряд Ук и прогноз ~к

1

0,8

0,6

0,4

0,2

0

На рис. 6 и 7 соответственно представлены оценки параметра модели (11) рк и частоты сезонной

составляющей ряда О к , полученные на основании алгоритма (19). Очевидно, что значение параметра в, описывающего сезонную составляющую ряда

у к, должно быть равно в « 0,6427. Однако оцениваемое значение в незначительно варьирует. Это 1,5

1

0,5 0

-0,5 -1

Рис.6. Оценка параметра Рк

-Х\ -

Рис.5. Средний квадрат ошибки прогнозирования

-2 1 к ,~ У Л

ек = YZ (~i - ~i) к i=1

76

РИ, 1998, № 3

связано с действующим шумом. Заметим, что в алгоритме (19) параметр Y был выбран равным 0,9 для усиления фильтрующих свойств алгоритма.

Далее рассмотрим процедуру идентификации параметров модели (3), (4), содержащей m гармоник неизвестных частот. Путём несложных преобразований можно привести (4) к виду

m-1

ук = Ё Р j+1(ук+j-m + Ук- j-m ) j=0

- У к-2m + £к = Pi2 yk - m + Рі(У к+1-m +

+ ук-1-m ) + в3(ук+2-m + ук-2-m ) + ••• + (2i)

+ Pm (ук-1 + ук-2m+1) - ук-2m + £к =

= рТу(к, m) - у ^2m + £к ,

(здесь Cj =--------), могут быть найдены путем

m j!(m-j)!

отыскания m корней степенного полинома от аргумента cos o .

Переходя к работе в реальном времени, записываем экспоненциально взвешенную процедуру стохастической аппроксимации, которая в данном случае имеет вид

Pn+1 = Pn + rN+1(.Pn+1 + уы-2m+1

< - PTNу(N +1 m))у(N +1 m),

rN+1 = Yn +|у(N +1,m)|2, 0 <y< 1.

При у = 1 она принимает форму алгоритма Гудвина-Рэмеджа-Кэйнеса, широко распространенного в адаптивных системах управления, а при у = 0 — алгоритма Уидроу-Хоффа, применяемого для обучения искусственных нейронных сетей.

При этом прогноз может быть записан в виде

у N+1 = Ріу(N + 1 m) - у N - 2m+1, (28)

а отфильтрованный сигнал —

УN+1,N+1 = РN+1у(N + 1, m) - ук-2m+1 • (29)

Определение оценок частот связано с необходимостью решения на каждом такте уравнения

где в = (Р1, Р2 ,•••, Pm )Т , У(k, m) = ^Ук-m ,

у к-m+1 + У к-m-1, У к-m+2 + У к-m-2 ,•••, У к-1 +

ук-2m+1)Т - (m X 1) векторы параметров и предыстории последовательности соответственно.

Вводя далее критерий идентификации

N ~Т 2

JI = Ё ((Ук + У к-2m ) -Р У(к, m)) (22)

к=2m+1

и используя стандартную процедуру наименьших квадратов, получаем вектор оценок Д. Далее можно

построить отфильтрованный сигнал (гармонический трецд)

УкN = РТ У(k, m) - У к-2m (23)

и прогноз

у N+1 = РТ У(N + 1 m) - У N - 2m +1 • (24)

Неизвестные частоты О j связаны с параметрами

в j соотношением

m-1

в1 + Ё в j+1 cos jo = cos mo (25)

и с учётом того, что

m /"*2 m-2 • 2

cos mo = cos o - Cm cos o sin o +

m-4 . 4 (26)

+ Cm cos o sin4 o + •••

^ m-1

P\, N+1 + Ё в j+1,N+1 cos jo N+1

j=1

= cos mo N+1,

(30)

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

Приведем пример, демонстрирующий работу по -лученного алгоритма. Ряд генерируется с помощью модели вида

ук = 0,93 sin 2п50к +1,34 sin 2п120к + £к •

Характеристики £ к следующие: математическое

ожидание равно 0 и <У ^ = 1,1. Ряд Ук и прогноз у к,

а также средний квадрат ошибки прогнозирования приведены соответственно на рис. 8 и 9.

4

3 2 1 0 -1 -2 -3

Рис.8. Ряд ук и его прогноз у к

РИ, 1998, № 3

77

2 1,5 1

0,5 0

На рис. 10 изображены оценки коэффициентов

в k и k , настраиваемые в соответствии с алгоритмом (27). Как и в предыдущем примере, в (27) Y равно 0 , 9.

- -Л /Кл. ^\/\ —

Рис.9. Средний квадрат ошибки прогнозирования ek

0,6 0,5 0,4 0,3 0,2

0,1 0

Рис. 10. Оценки коэффициентов в k и в2 k

Заметим, что из-за действия шума в спектре ряда на частоте / = 270 присутствует небольшой всплеск, что указывает на наличие ложной гармоники небольшой амплитуды.

В пользу адаптивного подхода к рассматриваемой задаче, кроме численной простоты, свидетельствует также следующее:

1. Возможность работы в условиях нестационарности, т. е. девиации частот. При этом следящие свойства алгоритма определяются варьированием параметра Y в (27).

2. Возможность работы в условиях “цветных”

шумов £k . При этом вектор y(k, m) расширяется

путём введения обновлений.

3. Защита от проблем, связанных с вырождением процедуры при больших N [12]. При этом объем обрабатываемой информации ограничивается в результате использования у < 1.

4. Возможность организации в реальном времени процедуры контроля над изменениями свойств обрабатываемой последовательности, т. е. раннего обнаружения разладок [8, 9, 16, 17].

5. Возможность распараллеливания вычислений и использования для обработки информации технологии искусственных нейронных сетей [18-20].

На рис. 11, 12 приведена схема искусственной нейронной сети, реализующей процесс прогнозирования и фильтрации тренд-сезонного временного ряда. Данная сеть состоит из двух независимых частей, первая

из которых (рис.11) обрабатывает полиномиальную компоненту ряда, а вторая (рис. 12) — сезонную. Входной слой сети на рис. 11 образован элементами

чистой задержки z -* 1, вследствие чего на первый скрытый слой, образованный сумматорами 2, параллельно подаются сигналы yk, Yk-1,• •, Yk_p . На выходе первого скрытого слоя формируются сигналы

разностей , V~k _!,•••> V~k - p+1, которые подаются на сумматоры второго скрытого слоя, вычисляющего вторые разности V 2 3 4 5~k, V 2~k-1,•••, V2~k_p+2. На выходе p -го скрытого слоя, образованного одним

сумматором, вычисляется сигнал p -й разности V p Yk .

Сигналы разностей подаются на усредняющие элементы а, реализующие рекуррентный алгоритм (8) и

вычисляющие средние значения yk, V~k,

V2~k, •••, VpYk . Сигналы среднего подаются на

релейные элементы, имеющие два состояния: 0, если соответствующее среднее равно нулю, и 1 в противном случае. Выходы реле подаются на входы логических

элементов Л, реализующих операцию типа

Ш = (x л у) л x и фактически фиксирующих наименьший порядок разности, не имеющей полиномиального тренда. Схема, приведенная на рис. 11, соответствует ряду с линейным трендом, у которого

V 2 ~k = 0 , что обнаружено логическим элементом,

на входы которого поданы сигналы x = 1, У = 0 . Контроль над выходами логических элементов позволяет обнаруживать в реальном времени изменения порядка полиномиального тренда. Кроме того, логические элементы управляют ключами К, которые открываются только при подаче на них сигнала, соответствующего единице. Таким образом, на входы

выходного сумматора 2 подаётся сигнал лишь той разности, которой соответствует единичный выход

логического элемента Л. В данном случае отпирается только ключ, соответствующий второй разности, в результате чего на выходе данной части сети появляется центрированный сигнал yk = V 2~k .

Вторая половина сети (рис. 12) обрабатывает сигнал Yk с исключенной полиномиальной компонентой и содержит во входном слое 2m элементов чистой задержки. Сумматорами первого скрытого слоя формируется вектор предыстории у(k, m) . Он подаётся на m -входовую адалину, обучение которой осуществляется с помощью алгоритма (27). На выходе адалины появляется отфильтрованная оценка уk, которая, будучи обработана цепочкой, состоящей из

p диграторов, даёт прогноз исходного ряда ~k • Дигрирование удобно также производить в соответствии с выражением

78

РИ, 1998, № 3

Рис. 11. Искусственная нейронная сеть для анализа полиномиального тренда

Рис. 12. Искусственная нейронная сеть для фильтрации тренд-сезонного временного ряда

РИ, 1998, № 3

79

уk = 1(-1)'Ciуk- -vpуk. (зі)

i=1

Контроль над возможными изменениями ряда y к осуществляется путём анализа поведения синаптических весов адалины в1, J32,..., вm и является

достаточно элементарной процедурой [9].

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

Литература: 1. Бокс Дж, Дженкинс Г. Анализ временных рядов. Прогноз и управление. М.: Мир, 1974. 406 с. 2. Чуев Ю.В., Михайлов Ю.Б., Кузьмин В.И. Прогнозирование количественных характеристик процессов. М.: Сов. радио, 1975. 400 с. 3. Кобринский Н.Е. Информационные фильтры в экономике. М.: Статистика, 1978. 287 с. 4. Бриллинджер Д. Временные ряды. Обработка данных и теория. М.: Мир, 1980. 536 с. 5. Кашьяп Р.Л., Рао А.Р. Построение динамических стохастических моделей по экспериментальным данным. М.: Наука, 1983. 384 с. 6. Льюис К.Д. Методы прогнозирования экономических показателей. М.: Финансы и статистика, 1986. 133 с. 7. Montgomery D.C., JohnsonL.A., GardinerJ.S. Forecasting and Time Series Analysis. N.Y.: McGraw-Hill. Inc., 1990. 384 p. 8. Isermann R. Fault diagnosis of machines via parameter estimation and knowledge processing. Tutorial paper // Automatica. 1993. 29. № 4. P.815-835. 9. Pouliezos A.D, Stavrakakis G.S. Real Time Fault Monitoring of Industrial Processes. Dordrecht: Kluwer Academic Publisher, 1994. 542 p. 10. Ивахненко А.Г., Мюллер Й.А. Самоорганизация прогнозирующих моделей. К.: Техніка, 1985; Берлин: ФЕБ Ферлаг Техник, 1984. 223 с. 11. Юрачковский Ю.П., Попков Н.В. Оценивание параметров в алгоритмах МГУА моделирования полигармонических процессов и полей // Автоматика. 1986. №6. С.9-16. 12. Shelekhova V. Yu. Harmonic algorithm GMDH for large data volume // SAMS. 1995. 20. P.117-126. 13. Бодянский Е.В., Плисс И.П, Соловьева ТВ. Многошаговые оптимальные упредители многомерных нестационарных стохастических процессов // Докл. АН УССР. 1986. Сер. А. №12. С.47-49. 14. Бодянский Е.В, Плисс И.П., Соловьева Т.В. Синтез квазипрямых адаптив-

ных регуляторов // Докл. АН УССР. 1987. Сер. А. №1. С.59-61. 15. Бодянский Е.В., Руденко О.Г. Адаптивные модели в системах управления техническими объектами. К.: УМК ВО, 1988. 212 с. 16. БодянскийЕ.В., Воробьёв С.А., Пл исс И.П. Адаптивная диагностика динамического объекта с периодическим выходным сигналом // Праці 3-ї Української конференції з автоматичного керування “Автоматика-96”. Т.1. Севастополь: СевГТУ, 1996. С.58-59. 17. Воробьёв С.А., Плисс И.П. Адаптивная диагностика динамического объекта с периодическим выходным сигналом // Деп. в УкрИНТЭИ 18.11.96. № 130. Харьков, 1996. 14 с. 18. Cichocki A., Unbehauen R. Neural Networks for Optimization and Signal Processing. Stuttgart: Teubner, 1993. 526 p. 19. Pham D.T., Liu X. Neural Networks for Identification, Prediction and Control. London: Springer-Verlag, 1995. 238p. 20. Scherer A. Neuronale Netze. Grundlagen und Anwendungen. Braunschweig/Wiesbaden: Friedr. Vieweg & Sohn Verlagsgesellschaft mbH, 1997. 249 s.

Поступила в редколлегию 23.07.98 Рецензент: д-р техн. наук Алексеев О.П.

Бодянский Евгений Владимирович, д-р техн. наук, профессор кафедры ТК ХТУРЭ. Научные интересы: теория адаптивных систем, искусственные нейронные сети, техническая диагностика, гармонический анализ. Хобби: фелинология, восточные учения, японская поэзия, история. Адрес: Украина, 310145, Харьков, пр. Ленина, 14. e-mail: [email protected]

Воробьев Сергей Анатольевич, канд. техн. наук, старший научный сотрудник ПНИЛАСУ ХТУРЭ. Научные интересы: искусственные нейронные сети, фильтрация и прогнозирование нестационарных процессов, фракталы и фрактальная размерность. Хобби: психология, иностранные языки, музыка. Адрес: Украина, 310726, Харьков, пр. Ленина, 14. e-mail: [email protected]

Костюк Ольга Васильевна, стажёр-исследователь кафедры САУ ХГПУ. Научные интересы: моделирование динамических систем, гармонический анализ, анализ нестационарных периодических процессов. Хобби: шитьё, горный туризм. Адрес: Украина, 310002, Харьков, ул. Фрунзе, 21.

Любчик Леонид Михайлович, д-р техн. наук, профессор кафедры САУ ХГПУ. Научные интересы: теория управления, адаптивные системы, анализ и прогнозирование случайных процессов, искусственные нейронные сети. Хобби: литература, искусство, живопись. Адрес: Украина, 310002, Харьков, ул. Фрунзе, 21. E-mail: [email protected]

80

РИ, 1998, № 3