Разработка методики интеграции формальных методов прогнозирования временных рядов и метода ассимиляции данных тема диссертации и автореферата по ВАК РФ 00.00.00, кандидат наук Тимошенкова Юлия Сергеевна

• Тимошенкова Юлия Сергеевна
• кандидат наук
• 2022

Оглавление диссертации кандидат наук Тимошенкова Юлия Сергеевна

Введение

Глава 1. Анализ предметной ситуации. Постановка исследования

1.1 Актуальность задачи прогнозирования

1.2 Основные термины и определения. Постановка задачи прогнозирования эволюции состояния динамической системы

1.3 Анализ формальных методов прогнозирования ВР, порожденных

ДС, описывающихся феноменологическими моделями

1.3.1 Авторегрессионные модели

1.3.2 Метод Singular Spectrum Analysis

1.3.3 Прогнозирование ВР на основе метода группового учета аргументов

1.3.4 Искусственные нейронные сети

1.3.5 Метод коррекции прогнозов при использовании формальных методах прогнозирования ВР

1.3.6 Итоги анализа формальных методов прогнозирования

1.4 Классический метод коррекции прогноза Data Assimilation

1.4.1 Фильтр Калмана

1.4.2 Ансамблевый фильтр Калмана

1.5 Постановка задач исследования

Глава 2. Разработка методики интеграции формальных методов

прогнозирования ВР и метода DA

2.1 Обоснование возможности интеграции формальных методов прогнозирования ВР и метода DA

2.2 Методика интеграции формальных методов прогнозирования ВР

и метода DА

2.3 Метрики, использованные для оценивания точности прогнозов значений ВР

2.4 Анализ точности прогнозирования СДУ Лоренца с помощью классического DA и предложенной методики

2.5 Выводы

Глава 3. Программный комплекс прогнозирования ВР

3.1 Описание программного комплекса

3.2 Описание ПБ № 1 «Формальные методы прогнозирования временных рядов»

3.3 Описание ПБ № 2 «Прогнозирование временных рядов с

помощью нейронной сети LSTM»

3.4 Описание ПБ № 3 «Прогнозирование временных рядов на основе интеграции формальных методов и метода ассимиляции данных»

3.5 Методы, входящие в состав вспомогательных модулей ПК ПВР

3.6 Результат использования программного комплекса ПВР

3.7 Выводы

Глава 4. Примеры использования программного комплекса

«Прогнозирование временных рядов»

4.1 Анализ точности прогнозирования набора данных пассажиропотока «Air Passengers»

4.1.1 Описание набора данных пассажиропотока «Air Passengers»

4.1.2 Анализ результатов прогнозирования ВР «Air Passengers»

4.2 Анализ точности прогнозирования набора данных солнечной активности

4.2.1 Описание набора данных числа солнечных пятен

4.2.2 Анализ результатов прогнозирования ВР солнечной активности

4.3 Анализ точности прогнозирования курса валют доллар/рубль и евро/рубль

4.3.1 Описание наборов данных курсов валют доллар/рубль и евро/рубль

4.3.2 Анализ результатов прогнозирования валютных курсов доллар/рубль, евро/рубль

4.4 Выводы

Заключение

Список сокращений и условных обозначений

Список литературы

Приложение А. Копии актов о внедрении

Приложение Б. Свидетельства о регистрации программы для ЭВМ

Введение

Введение диссертации (часть автореферата) на тему «Разработка методики интеграции формальных методов прогнозирования временных рядов и метода ассимиляции данных»

Актуальность темы.

Под термином «прогноз», традиционно, понимают научно обоснованное суждение о возможных состояниях выбранного объекта (динамической системы (ДС)) в будущем и/или его (ее) траекториях в пространстве состояний, а также сроках осуществления данных состояний. Метод прогнозирования, помимо общеизвестного ежедневного прогноза погоды, применяется, например: в автоматических системах управления производственными процессами [1]; для прогнозирования геофизических и метеорологических явлений на Земле [2; 3]; в геологии и геофизике [4; 5]; в медицине [6], [7—9]; в политологии [10]; в макроэкономике и государственном управлении [11—14], в экономике [15], в управлении финансами [13; 15]

Таким образом, задача прогнозирования будущих состояний ДС на основе информации о ее предыдущих состояниях является актуальной в различных областях человеческой деятельности [16].

Во многих случаях задачу прогнозирования изменения состояния ДС удается свести к задаче прогнозирования упорядоченных по времени последовательностей, составленных из значений параметров процессов X (Ь), порожденных данной ДС, измеренных в узлах равномерной временной сетки Ь = гАЬ, г = 0, N — 1, где АЬ = сош^ шаг временной сетки. Данные последовательности называют временным рядом (ВР). При этом различают следующие типы ДС:

- ДС, описываемые известными математическими моделями, которые построены на основе физических законов, определяющих ее эволюцию во времени;

- ДС, описываемые с помощью феноменологических моделей порожденных ими ВР - математических моделей ВР, построенных на основе их анализа.

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

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

Степень разработанности темы.

Начало разработки формальных методов прогнозирования ВР датируется 1910 г, когда было предложено использовать для прогнозирования ВР модель скользящего среднего (MA). В 1920 г в работах G. Yule и J.C. Walker впервые была предложена авторегрессионная модель (AR) ВР. Далее в 1970 г. G. Box и H. Jenkins объединили AR- и MA-модели ВР в единую модель, получившей название авторегрессии-скользящего среднего (ARMA). Затем в 1976 г. на основе ARMA-модели G. Box и H. Jenkins была разработана модель авторегрессии-про-интегрированного скользящего среднего (ARIMA).

В 1989 г R. Vautard, M. Ghil предложили метод сингулярного спектрального анализа (SSA), относящийся к классу непараметрических методов анализа и прогнозирования ВР, основанный на декомпозиции исследуемого ВР на некоторый аддитивный набор более простых ВР, называемых главными компонентами (ГК), которые могут быть интерпретированы с точки зрения представлений о динамике ДС, породившей данный ВР (например, тренд, гармонические составляющие, шум и т.д.). В основании метода SSA лежит преобразование одномерного ВР в многомерный ВР с помощью однопараметрической сдвиговой процедуры (от-

сюда в русскоязычной литературе данный метод получил название «Гусеница») и далее исследование полученной многомерной траектории с помощью сингулярного разложения и восстановления (аппроксимации) ВР по выбранным ГК. Отметим, что данный метод не требует стационарности анализируемого ВР, в отличие, например, от Фурье-анализа и параметрических методов спектрального оценивания. Исследованием различных теоретических и прикладных аспектов метода SSA также занимались Н.Э. Голяндина, А.А. Жиглявский, С.В. Поршнев и др.

В 1968 г А.Г. Ивахненко был предложен формальный метод прогнозирования ВР, получивший название метод группового учета аргументов (GMDH), основанный на гипотезе о том, что значение xN ВР xi зависит от B его предыдущих значений:

XN = F (XN-i,XN-2, • • • ,XN-B) . (1)

Здесь в качестве функции F (...) он предложил использовать полином Кол-могорова-Габора:

m mm

XN = ao + J2 aiXN-i + J2 ai,JXN-iXN-j + i=l i=l j=l

m m m

ai,j k XN-iXN-j XN-k,

i=1 j=1 k=1

m ^ B, коэффициенты которого вычисляются с использованием метода наименьших квадратов (МНК), который, однако, оказывается весьма затратным с вычислительной точки зрения.

Далее в 1982 г им же был предложен самоорганизующийся GMDH, в котором используются полиномы Колмогорова-Габора второго порядка, аргументами которых являются все возможные сочетания пар входных переменных {xN-1, xN-2,...,xN-B}. В результате строится множество регрессионных моделей ВР первого уровня. Из них выбираются те регрессионные модели, дисперсия ошибок остатков которых оказывается меньше выбранного значения. Затем строятся значения, в которых в качестве входных переменных полиномов используются регрессионные модели второго уровня полиномов первого уровня т.д. Описанная выше процедура построения полиномов продолжается до тех пор пока погрешность полиномов данного уровня не оказывается меньше погрешности полинома предыдущего уровня. Данный алгоритм оказывается менее затратным с вычислительной точки зрения, нежели классический GMDH.

Один из активно используемых в настоящее время формальных методов прогнозирования ВР основан на использовании искусственных нейронных сетей (1943 r. W. McCulloch и W. Pitts), выполняющих роль функции F (...) в (1.13). При этом оказывается, что наилучшее качество прогнозов, в рассматриваемой задаче, обеспечивают рекуррентные нейронные сети (RNN) с долгой краткосрочной памятью (LSTM) (1997 г. J. Hochreiter и J. Schmidhuber).

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

Свободным от отмеченного выше недостатка формальных методов прогнозирования ВР оказывается метод ассимиляции данных (Data Assimilation, DA), предложенный в работах M. Talagrand и J. Courier. Однако, обязательным условием использования DA является знание математической модели ДС, породившей данный ВР, которое, как оказывается, выполняется на практике далеко не всегда.

Целью работы является разработка методики интеграции формальных методов прогнозирования ВР, порожденных ДС, описываемыми феноменологическими моделями, и метода Data Assimilation.

Для достижения поставленной цели потребовалось решить следующие задачи:

1. Анализ существующих методов прогнозирования и коррекции ВР.

2. Разработка методики интеграции формальных методов прогнозирования ВР в DA.

3. Разработка программной реализации методики интеграции формальных методов прогнозирования ВР в DA.

4. Экспериментальная апробация разработанной методики интеграции формальных методов прогнозирования ВР в DA на примере ВР, порожденных ДС, описываемыми феноменологическими моделями.

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

Предмет исследования: методика интеграции формальных методов прогнозирования ВР, порожденных ДС, описываемыми феноменологическими моделями, и метода DA.

Научная новизна: диссертационного исследования заключается в:

1. Обосновании возможности интеграции формальных методов прогнозирования ВР, порожденных ДС, описываемыми феноменологическими моделями, и метода DA (соответствует п. 11 паспорта специальности 05.13.01 Системный анализ, управление и обработка информации).

2. Разработке методики интеграции формальных методов прогнозирования ВР в метод DA и ее программной реализации (соответствует п. 5, 11, 12 паспорта специальности 05.13.01 Системный анализ, управление и обработка информации).

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

Теоретическая и практическая значимость диссертационного исследования заключается в разработке методики интеграции формальных методов прогнозирования ВР, порожденных ДС, описываемыми феноменологическими моделями, и метода DA и программного комплекса, обеспечивающего использование данной методики для прогнозирования реальных ВР, в которой реализованы:

1. Формальные методы прогнозирования ВР: ARIMA, SSA, GMDH, LSTM-

net.

2. Классический метод DA.

3. Методика интеграции формальных методов прогнозирования ВР и метода Data Assimilation.

Работа выполнена в рамках проекта № 20-37-90006/20 от 21.08.2020, поддержанного грантом РФФИ «Аспиранты» 2020 г.

Методология и методы исследования. В работе использованы методы анализа и прогнозирования ВР, методы ассимиляции данных, методы математической статистики, вычислительной математики и системного анализа.

Основные положения, выносимые на защиту:

1. Разработана научно обоснованная методика, обеспечивающая интеграцию формальных методов прогнозирования ВР и метода DA.

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

3. Результаты прогнозирования ВР, представляющих собой ВР, порожденные ДС, описываемыми феноменологическими моделями, подтверждают работоспособность разработанной методики интеграции формальных методов прогнозирования ВР и метода DA и ее программной реализации, использованы: ВР «Air Passengers», составленный из данных о числе пассажиров, перевезенных авиакомпаниями США в период с января 1949 по декабрь 1960; ВР, составленного из среднемесячных значений чисел Вольфа (числа наблюдаемых солнечных пятен), сглаженных по 13 месяцам, являющихся количественными показателями активности Солнца с августа 1913 по сентябрь 2021 гг.; ВР, составленного из значений переводных курсов доллар/рубль, евро/рубль в период с ноября 2021 г. по апрель 2022 г.; представляющие собой ВР, порожденные ДС, описываемыми феноменологическими моделями, подтверждают работоспособность разработанной методики интеграции формальных методов прогнозирования ВР и метода DA и ее программной реализации.

Достоверность полученных результатов подтверждается обоснованным применением методов прогнозирования ВР, результатами тестовых испытаний разработанного программного комплекса «Прогнозирование временных рядов», а также сходимостью результатов прогнозирования ВР с помощью разработанной методики интеграции с соответствующими результатами их наблюдения.

Внедрение результатов диссертационного исследования. Результаты диссертационного исследования внедрены в ООО «Эйрбэйс» (акт об использовании результатов от 05.03.2022), а также используются в федеральном государственном автономном образовательном учреждении высшего образования «Уральский федеральный университет им. первого президента России Б.Н. Ельцина» при подготовке бакалавров и специалистов по укрупненной группе специальностей 10.00.00 Информационная безопасность (акт об использовании результатов от 21.07.2022).

Апробация работы. Результаты работы докладывались и обсуждались на научно-практических конференциях: Международная конференция студентов, аспирантов и молодых ученых «Информационные технологии, телекоммуникации и системы управления», 2017, Екатеринбург, Россия; XI International scientific and technical conference «Applied Mechanics and Dynamics Systems», 14-16 ноября 2017, Омск, Россия; 1st International Conference on Physics, Mathematics and Statistics (ICPMS 2018), 12-14 May 2018, Shanghai, China; 2018 International Conference on Applied Mathematics and Computational Science (ICAMCS.NET 2018),

6-8 october 2018, Budapest, Hungary; International Conference on Numerical Analysis and Applied Mathematics 2018 (ICNAAM 2018), 13-18 September 2018, Rhodes, Greece; 4th International Workshop on Radio Electronics and Information Technologies (REIT-Autumn 2018), Екатеринбург, Россия; Международная конференция 2019 Ural Symposium on Biomedical Engineering, Radioelectronics and Information Technology (USBEREIT), 25-26 апреля 2019, Екатеринбург, Россия; Международная конференция 2020 Ural Symposium on Biomedical Engineering, Radioelectronics and Information Technology (USBEREIT), 14-15 мая 2020, Екатеринбург, Россия; Международная конференция 2021 Ural Symposium on Biomedical Engineering, Radioelectronics and Information Technology (USBEREIT), 13-14 мая 2021, Екатеринбург, Россия.

Личный вклад.

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

Публикации. По теме диссертации опубликовано 11 научных работ, в том числе 9 научных статей в изданиях, определенных ВАК РФ и Аттестационным советом УрФУ, из них 8 статей в изданиях, индексируемых в международных цитатно-аналитических базах Scopus и Web of Science. Зарегистрированы 3 программы для ЭВМ.

Объем и структура работы. Диссертация состоит из введения, 4 глав, заключения и 2 приложений. Полный объём диссертации составляет 162 страницы, включая 72 рисунка и 6 таблиц. Список литературы содержит 105 наименований.

Глава 1. Анализ предметной ситуации. Постановка исследования 1.1 Актуальность задачи прогнозирования

Прогнозом (прогноз - от греч. проу^ст^ «предвидение, предсказание») будем называть, следуя [17], научно обоснованное суждение о возможных состояниях выбранного объекта (системы) в будущем и/или его (ее) траекториях движения в пространстве состояний, а также сроках осуществления данных состояний.

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

1. В автоматических системах управления (АСУ) производственными процессами (ПП), в которых на основе сравнения спрогнозированных и реально наблюдаемых параметров соответствующих ПП вырабатывается то или иное воздействие на технического устройство, управляющее данным ПП [1; 18].

2. При прогнозировании геофизических и метеорологических явлений на Земле (в том числе, погоды), обусловленных влиянием космических и гелиофи-зических (солнечная активность) факторов [2; 3].

3. В геологии и геофизике для оценивания запасов ископаемых ресурсов и горно-геологических условий бурения [4; 5].

4. В медицине для прогнозирования: распространения заболеваний [6] и эпидемий [7], продолжительности жизни пациента, когда известны его наследственные заболевания [8], уровни кислорода и сахара в крови, а также корректировки проводимого лечения на основе анализа отклонений спрогнозированных значений контролируемых показателей от их фактических значений [9].

5. В политологии для прогнозирования развития текущей социально-политической ситуации [10].

6. В макроэкономике и государственном управлении при стратегическом планировании, в том числе для прогноза объемов внутреннего валового продукта (ВВП) [11], динамики роста населения [12], инфляции, изменения количества рабочих мест, объемов производства, потребления товаров и ресурсов [13], а также обоснованного формирования налогово-бюджетной, кредитной политик,

планирования расходов на здравоохранение, пенсионные программы и программы борьбы с бедностью [14].

7. В экономике при планировании производства и его размещения, управления запасами, определения ассортимента продукции [15].

8. В управлении финансами, в том числе для прогнозирования динамики доходности тех или иных финансовых операций и стоимости ценных бумаг [13;

15].

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

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

В 1929 г Министерство труда США своем новогоднем прогнозе объявило, что 1930 г будет хорошим в плане трудоустройства граждан, однако, фактическая ситуация оказалась прямо противоположной - 1930 г стал годом мирового экономического кризиса, который привел экономику к стагнации и, как следствие, сокращению числа рабочих мест [20].

В 1756 г. был дан прогноз о том, что численность населения Земли останется неизменной [20]. Однако, в действительности оказалось, что численность населения планеты Земля, несмотря на войны, продолжает увеличиваться. При этом, по мнению некоторых экспертов (см, например, [21]) темпы роста численности населения с течением времени непрерывно увеличиваются.

В 1899 г. руководитель американского патентного бюро Чарльз Дьюел предложил закрыть патентные бюро, так как согласно его прогнозу все изобретения уже сделаны [20].

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

1. Методы, основанные на использовании экспертных оценок (класс № 1).

2. Методы, основанные на использовании математических моделей, описывающих динамику исследуемых объектов (систем) (класс № 2).

3. Методы, основанные на использовании феноменологических моделей, описывающих динамику исследуемых объектов (систем) (класс № 3).

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

При этом принимая во внимание, что методы прогнозирования на основе экспертных оценок ранее рассматривались в работах многих исследователей (см., например, [22; 23] и др.), объектом исследования выбраны методы прогнозирования относящиеся к классам № 2, 3.

1.2 Основные термины и определения. Постановка задачи прогнозирования эволюции состояния динамической системы

Введем основные термины и определения, используемые далее в нашем исследовании.

Определение 1. Динамической системой (ДС) будем назвать физический, биологический, экономический и др. объект, состояние которого изменяется во времени.

При этом будем различать следующие типы ДС:

1. ДС первого типа - ДС, описываемые известными математическими моделями (например, механические колебания, конвекция воздуха; колебания уровня мирового океана; температура поверхности моря и т.д.).

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

Определение 2. Состоянием ДС в момент времени Ь будем называть некоторую совокупность показателей, характеризующих ДС в данный момент времени:

х(г) = {хг(ь),х2(ь), ...,хп(г)}.

Определение 3. Пространством состояний ДС будем называть множество состояний ДС в известные моменты времени:

[X(ьо),х(и), ...,х(¿ж-1)}.

Определение 4. Проекцию пространства состояний ДС на выбранную координатную ось будем называть процессом .

Определение 5. Временным рядом (ВР) Х{, г = — 1 будем называть некоторую упорядоченную во времени выборку мгновенных значений процесса Х^), порожденного той или иной ДС:

{ж0 = х (г0),х1 = х (¿1),..., хм—1 = х (Ьм-1)},

где г0 < г1 < ... < гм—1.

Так как ВР является одномерной проекцией пространства состояний ДС на выбранную координатную ось, его свойства определяются законами эволюции данной ДС.

Определение 6. Равномерным ВР будем называть ВР, составленный из мгновенных значений процесса х(Ь), измеренных в узлах равномерной временной сетки = гАЬ, г = 0, N — 1, где АЬ =const- шаг временной сетки. С математической точки зрения равномерный ВР представляет собой некоторую упорядоченную последовательность действительных чисел.

Определение 7. Неравномерными ВР будем называть ВР, составленный из мгновенных значений процесса х (Ь), измеренных в узлах неравномерной временной сетки Ь0 < Ь1 < ... < —1 такой, что существуют такие точки, для которых выполняется условие Ьт+1 — Ьт = Ьт+2 — Ьт+1, т = 0,N — 3. С математической точки зрения неравномерный ВР представляет собой таблицу (Ь^ х^.

Очевидно, что неравномерный ВР может быть преобразован в равномерный ВР с помощью известных методов интерполяции функций, заданных таблично [24].

Год

Рисунок 1.1 — Временной ряд среднегодовых значений температуры в районе океанической станции погоды «INDIA» на глубине 200 м [25]

Пример равномерного ВР, содержащего реальные наблюдения среднегодовых значений температуры в районе океанической станции погоды «INDIA» [25], представлен на рисунке 1.1.

Отметим, что многомерные процессы X{t) = [Xi(t) ,Х2(t), ..порождают многомерные ВР (с математической точки зрения - матрицы):

Хи (to) Хи (¿i)

Xl( tN-i) Х2( tN-i) ••• ХИ (tN-l)

Для данного типа ВР задача анализа и прогнозирования, решается для каждого из столбцов матрицы X.

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

Задача прогнозирования ВР, порожденных ДС, описываемой с помощью известной математической модели (Задача № 1), имеет следующую формулировку.

X=

Xi(to)

Xi(ti)

X2(to)

X2(ti)

11,5

11

Р 10,5

ей

л

ь

10

л

<и'

с

М 9,5

н '

9

8,5

1960 1965 1970 1975 1980 1985 1990 1995 2000 2005 2010 2015

Год

Рисунок 1.2 — К постановке задачи прогнозирования ВР

Задача № 1. Известна математическая модель, описывающая эволюцию ДС в дискретном времени:

хк+1 = м (хк ,гк) + ъик, (1.1)

где X - вектор состояния ДС, размерностью Ь х 1, Ь = length(X); М - оператор или функция перехода, определяющая эволюцию ДС во времени; Хк - вектор состояния ДС в момент времени Ьк, к = 0,1,...,Ж — 2; Хк+х - вектор состояния в момент времени Ьк+\; ъхк - случайный вектор ошибок математической модели, у которого неизвестны точные значения его координат, но известны их статистические характеристики. Соответственно, в каждый момент времени Ьк, к = 0,1,...,N — 1 известны значения проекции вектора состояния X на выбранную координатную ось - ВР хк.

Требуется на основе имеющейся информации получить статистически устойчивые оценки значений ВР х^т, т = 1,2,..., соответствующие моментам времени Ь^—\+т.

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

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

1. Генерация согласно (1.1) ансамбля значений векторов состояния ДС {х}п, п = 0,1,...N — 1, длиной Ь = N — 1 + N5, N - число прогнозируемых точек (дальность прогноза) и вычисление далее их проекций на выбранную координатную ось - ансамбля ВР {Хг}п, г = 0,1,N — 1 + N, п = 0,1,...N — 1. (При вычислении указанных ансамблей в качестве шума математической модели ъик используются независимые случайные числа, генерируемые в соответствии с известными законами распределения.)

2. Формирование из ансамбля ВР {х}п, г = 0,1,N — 1 + N, п = 0,1,...,Да — 1 нового ансамбля |Х^т, т = 0,1,..М — 1, составленного из таких ВР, что

. N —1

^Е Х — Ш2 < а2,

¿=0

где А - выбранное значение погрешности.

3. Вычисление прогнозируемых значений ВР хт, т = 1,2,..., соответствующие моментам времени ЬМ—1+т, как среднего по ансамблю ., ] = 0,1,..М:

х

т

1 М—1 3 =0

В наиболее общем виде задача прогнозирования ВР, порожденного ДС, описываемой с помощью феноменологической модели (Задача № 2), имеет следующую формулировку.

Задача № 2. Известны значения ВР х^, г = 0^ — 1, порожденного ДС, описываемой с помощью феноменологической модели. Требуется на основе имеющейся информации вычислить прогнозируемые значения ВР х^т, т = 1,2,..., соответствующие моментам времени Ьт.

Список литературы диссертационного исследования кандидат наук Тимошенкова Юлия Сергеевна, 2022 год

тр - -

N

к=1

И • а = т, где И - автокорреляционная матрица:

1 п т2 • • • тр-1

п 1 т1 • • • т'р-2

И = т'2 т1 1 • • • тр-3

_ тр-1 т'р-2 тр-3 • • • 1

а1

«2

а = ,

ар

1

т\

т =

Тр-1

При использовании для прогнозирования AR(p)-модели ее порядок р совпадает с базой прогноза В, р = В.

Схема алгоритма прогнозирования ВР на основе AR-модели представлена на рисунке 1.6, из которого видно, что данный алгоритм реализуется выполнением следующей последовательности действий:

1. Проверить ВР на стационарность, используя для этого, например, KPSS-тест [26].

2. Если ВР не стационарный, то использовать известные методы приведения нестационарного ВР к стационарному ВР (например, вычитание трендовой составляющей, сезонной компоненты, логарифмированием ВР) и перейти в п.2, иначе перейти в п.3:

Рисунок 1.6 — Схема алгоритма прогнозирования с помощью AR-модели

3. Провести анализ ВР и выбрать порядок р (как правило, достаточно р ^

2).

4. Рассчитать параметры модели а1,... ,ар. Обычно для оценки параметров авторегрессионного процесса с использованием заданного временного ряда используются уравнения Юла—Уокера [27].

5. Вычислить прогнозируемое значение ВР по формуле:

р

Хм = с + ег + агХм_г.

¿=1

6. N = N + 1.

7. перейти в п.5.

Модель скользящего среднего

В модели скользящего среднего (Moving average, MA) q-го порядка, обозначаемой MA(q), связь между xN и xN-1,xN-2,..., xN-q описывается выражением:

q

XN = М + iN + в6N-i, (1-7)

i=1

где q - порядок модели, вь • • • в - параметры модели, ц - среднее значение ВР, ti - собой отсчеты белого шума (нормального случайного процесса с нулевым средним и постоянной дисперсией) [28]. Параметры модели MA(q) (см. (1.7)) находятся с использованием метода максимального правдоподобия (ММП) [28].

Отметим, что сегодня функции, возвращающие значения параметров как MA(q), так и других авторегрессионных моделей являются неотъемлемой составляющей каждого из современных математических пакетов и известных программных библиотеках языков программирования высокого уровня, в том числе и языка Python 3, поэтому здесь и далее не приводятся известные формулы, по которым рассчитывают значения оцениваемых параметров авторегрессионных моделей.

Схема алгоритма прогнозирования на основе модели MA(q) представлена на рисунке 1.7.

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

1. Проверка ВР на стационарность, например, с помощью KPSS-теста [26].

2. Если ВР не стационарный, то ВР приводится к стационарному ВР (например, с помощью удаления трендовой составляющей, сезонной компонент, логарифмирования ВР, перехода от исходного ВР к ВР, содержащего значения разностей d-го порядка), иначе перейти кп.3.

3. Оценка порядка q MA^-модели.

4. Вычисление оценок параметры MA(q)-модели вь • • • в с помощью метода наименьших квадратов или ММП.

5. Вычисление прогнозируемых значений ВР согласно (1.7).

6. N = N + 1.

7. перейти в п.5.

Рисунок 1.7 — Схема алгоритма прогноза с помощью MA

С вычислительной точки зрения оценка параметров МЛ(д)-модели ВР является весьма затратной процедурой, поэтому для прогнозирования ВР, порожденных реальными ДС, МА(д)-модели ВР пользуется крайне редко.

Модель авторегрессии проинтегрированного скользящего среднего

Модель авторегрессии проинтегрированного скользящего среднего (ARIMA - Autoregressive Integrated Moving Average) используется для описания ВР xi, значения которых вычисляются по следующей формуле [28]:

р я

хм = С + ем + - + ^ вз ем -3, (1.8)

¿=1 з=1

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

Также модель ARIMA, используя лаговый оператор (оператор смещения), можно представить Ь : Ьхм = хМ1 в виде:

р я

(1 - Ь)Х = с аЬг)(1 - Ь)Х + (1 + ^ взЬз)ем-1.

¿=1 з=1

Оценки значений параметров модели ARIMA вычисляются выполнением следующей последовательности действий:

1. Оценка порядка ё на основе проверки стационарности ВР: если исходный ВР стационарный - ё = 0, иначе вычисляются первые разности исходного ВР (Ьх,) и далее проверяют ВР Ьх, на стационарность: если ВР Ьх, нестационарный, вычисляют вторые разности ВР х, и т.д. (Для оценки порядка ё также можно использовать метод единичных корней [16].)

2. Вычисление AR-параметров а1, ..., ар.

3. Вычисление МА-параметров в1, • • •, вя.

Схема алгоритма прогнозирования ВР х, с помощью метода ARIMA представлена на рисунке 1.8, из которого видно, что он реализуется выполнением следующей последовательности действий.

1. Извлечение значения ВР х, из хранилища данных.

2. Проверка ВР х, на стационарность, например с помощью KPSS-теста [26].

3. Если ВР х, стационарный, то установить ё = 0 и перейти к п. 4, иначе установить ё =1 и перейти от ВР х, к ВР Ьх,, составленному из первых разностей, аналогично, проверить его на стационарность: если ВР Ьх, стационарный ВР перейти к п 4, иначе установить ё =1 перейти от ВР Ьх, к ВР Ь2х, и т.д.

4. Оценить на основе ВР порядка ё параметры р, д ARIMA-модели.

5. Вычислить оценки параметров ARIMA-моделей а1,... ,ар, в1, • • • ,вя.

Рисунок 1.8 — Схема алгоритма прогноза с помощью ARIMA

6. Вычислить прогнозируемое значение ВР хм согласно (1.8).

7. Вычислить будущее значение ВР хм согласно (1.8).

8. N = N + 1.

9. перейти в п.6.

1.3.2 Метод Singular Spectrum Analysis

Метод сингулярного спектрального анализа (англ. - Singular Spectrum Analysis, SSA) [5; 29; 30] является непараметрическим методом анализа и прогнозирования ВР. Он основан на декомпозиции исследуемого ВР на некоторый аддитивный набор более простых ВР, называемых главными компонентами (ГК), которые могут быть интерпретированы с точки зрения представлений о динамике ДС, породившей данный ВР (например, тренд, гармонические составляющие, шум и т.д.). Данный метод основан на преобразовании одномерного ВР в многомерный ВР с помощью однопараметрической сдвиговой процедуры (отсюда в русскоязычной литературе данный метод получил название «Гусеница») и далее исследовании полученной многомерной траектории с помощью сингулярного разложения и восстановления (аппроксимации) ВР по выбранным ГК. Отметим, что метод SSA не требует стационарности анализируемого ВР, в отличие, например, от Фурье-анализа и параметрических методов спектрального оценивания.

Приведем, следуя [29], краткое описание базового алгоритма SSA на примере вещественнозначного ВР XN = (x0, x2,..., xN _ 1) длинной N. При этом будем считать, что N > 2. Предположим, что ряд XN ненулевой, то есть существует как минимум одно i, такое, что Xi = 0. Базовый алгоритм состоит из двух взаимодополняющих друг друга этапов - этапа разложения и этапа восстановления.

Этап 1 - Разложение:

Шаг 1. Вложение - На данном шаге формируется траекторная матрица. Пусть L - некоторое целое число, 0 < L < N. Процедура вложения образует K = N _ L + 1 векторов вложения:

Xi = (Xi—1, . . . ,Xi+L_2) , 1 < i < K,

имеющих размерность L. Если нужно подчеркнуть размерность векторов Xi, то их называют векторами L-вложения.

Матрица X, составленная из векторов L-вложения в качестве столбцов

X = [Xi : ... Xk ],

называется траекторной матрицей.

Другими словами, траекторная матрица - это матрица:

x\ x2 ... xK

X= Xj = (x2 X3 7x*+1 )■ (19)

XL xl+1 ... xT

Очевидно, что Xj = xi+j-_2, поэтому матрица (X) имеет одинаковые элементы на диагоналях i + j =const, то есть является ганкелевой. Существует взаимно-однозначное соответствие между ганкелевыми матрицами размерности L х K и ВР длины N х L х K х 1.

Шаг 2 - Разложение по сингулярным числам (Singular Value Decomposition, SVD).

На данном шаге выполняется SVD траекторной матрицы.

Пусть S = xxt. Обозначим собственные числа матрицы S, взятые в неубывающем порядке, Аь ..., XL, А1 ^ ... ^ XL и U1, ..., UL - соответствующую ортонормированную систему собственных векторов матрицы S.

Пусть d = rankx = max {i : Аi > 0}. Если обозначить Vi = xtUi/y/Ai, i = 1, ..., d, то результат сингулярного разложения матрицы X может быть записан как:

X = Xi + ... + xd, (1.10)

где xi = V\UVtT. Каждая из матриц xi имеет ранг 1, поэтому их можно назвать элементарными матрицами. Набор {^/\i, Ui, VT} называют ¿-ой собственной тройкой.

Этап 2 - Восстановление (группировка и усреднение):

Шаг 3 - Группировка. - На основе разложения (1.10) процедура группировки делить все множество индексов 1, ..., d на m непересекающихся подмножеств I1, ..., Im. Пусть I = i1, ..., ip, тогда результирующая матрица Xj, соответствующая группе I определяется как:

XJ = Xi1 + ... + Xip.

Группировкой называется процедура выбора множеств I = I1, ..., Im.

Шаг 4 - Диагональное усреднение. - На данном шаге каждая матрица размером L на K сгруппированного разложения переводится в новый ряд длины N.

Пусть Z = xj = xjl + ... + xim некоторая матрица размером L на K с элементами zij, где 1 ^ i ^ L, 1 ^ j ^ K. Предположим L* = min (L,K),

К * = тах (Ь,К) и N = К + Ь — 1. Пусть ^ = , если Ь < К, и ^ = иначе.

Диагональное усреднение переводит матрицу Ъ во ВР д0, ..., дт—1 по формуле:

{

1 2* 0 < к < Ь* — 1

к+1 ш=1 2т,к—т+2 0 < к < Ь 1,

= ' £ Ет=1 к—т+2 Ь* — 1 < к < К*, (111)

1 VМ—К ф+1 2 * К * < к < N

М-к 2^т=к—К*+2 2т,к—т+2 К < к < N'

gfc = \

где L* = min(L,K), k* = max(L,K) и X = L + K — 1.

Выражение (1.11) соответствует усреднению элементов матрицы вдоль диагоналей, на которых i + j = k + 2 : выбор k = 0 дает g0 = уц, для k = 1 получаем

gi = (Ш + У2)/2 и т.д.

Заметим, что, если матрица Y является траекторной матрицей некоторого ВР (h0,... ,hN) (другими словами, если матрица Y является ганкелевой), то gi = hi для всех i. Применяя диагональное усреднение (1.11) к результирующим матрицам, получаем ВР X = (Х0, ... ,xN—1). Следовательно, исходный ряд X = (x0,... ,xN — 1) раскладывается на сумму m рядов:

m

x = £ xk. k=i

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

Метод SSA как средство прогнозирования

Метод SSA также может быть использован для прогнозирования ВР. Сегодня известны следующие алгоритмы прогнозирования ВР, использующие данный метод:

- Recurrent SSA (SSA-R);

- Stochastic SSA (SSA-S).

Алгоритм SSA-R реализуется выполнением следующей последовательности действий:

1. Извлечение значения ВР г = 0,Ы — 1 из хранилища данных.

2. Задание числа столбцов траекторной матрицы Ь.

3. Формирование траекторной матрицы X ВР Х{, г = 0,Ы — 1, согласно (1.9).

4. Вычисление матрицы S = ххт.

5. Вычисление собственных чисел Х и собственных векторов Щ, г = 1,Ь матрицы S.

6. Вычисление значений векторов У = ХТи^у/Х^, г = 1,Ь. (Результат -набор собственных троек {\Г1, Щ, У}), сингулярного разложения траекторной матрицы X.)

7. Выбор номеров собственных троек, используемых для прогнозирования значений ВР х^. (Результат-вектор I = 11,... ,1т, содержащий номера выбранных собственных троек.)

8. Вычисление матриц Х1х,..., Х1т,

10. Диагональное усреднение матрицы Ъ согласно (1.11) (Результат - ВР дк,

11. Задание длины прогноза М.

12. Выбор из матрицы и1к, составленной из собственных векторов У1к, последнего столбца матрицы и^.

13. Вычисление нормы V вектора и^м.

14. Вычисление весовых коэффициентов Я.

Х1к = у/ХКи1к УТ,

к = 1,... ,т.

9. Вычисление матрицы Ъ:

т

¿=1

к = 1,Х).

(1.12)

15. Вычисление спрогнозированного значения ВР хм:

Ь—2

хм = Е дм-3—1-

3=0

16. N = N + 1.

17. Повторить п.3-16 для ВР х1, ... ,ХМ.

Алгоритм SSA-S реализуется выполнением следующей последовательностью действий:

1. Извлечение значения ВР XI, г = 0^ — 1 из хранилища данных.

2. SSA-анализ ВР с целью выбора параметра Ь и группируемых собственных троек Х1к, и1к, Ухк сингулярного разложения ковариационной матрицы анализируемого ВР.

3. 3 = 0. _

4. Добавление к ВР XI, г = 0^ — 1 отсчета, являющегося случайным числом, сгенерированным в соответствии с равномерным законом распределения на интервале [тт(х^); тах(х^)] (результат - ВР х^ г = 0,N).

5. Декомпозиция ВР, полученного в п. 4 ВР х i, г —

6. Группировка вычисленных собственных троек ВР XI, г = 0^, номера которых выбраны в п. 2.

7. Восстановление ВР х^ г = 0^.

8. Хм = хм.

9. Если 3 = 0, то выполнить п. 11, иначе перейти в п. 10.

10. Если \хМ} — хМ 1}| < е, е - выбранное значение погрешности, то перейти к п. 13, иначе перейти к п. 11.

11. 3 = 3 + 1

12. Повторение пп. 5-8.

13. Выбор последнего отсчета как значения прогноза.

14. N = N + 1.

15. Переход к п. 4.

Схемы данных алгоритмов представлены на рисунках 1.9 и 1.10, соответственно.

Данный метод не требует строгой стационарности ВР и успешно применяется для прогнозирования ВР различной природы [32—35]: экономические ВР, промышленные и социальные ВР, в том числе [29] для прогноза ВР, порожденных ДС, описываемых с помощью феноменологических моделей, например, объем пассажироперевозок международным аэропортом.

Pисунок 1.9 — Схема алгоритма прогноза с помощью SSA-R

Рисунок 1.10 — Схема алгоритма прогноза с помощью SSA-S

Главным преимуществом метода SSA-S в сравнении с методом SSA-R состоит в том, что первом методе выбор номеров группируемых собственных троек осуществляется на первом шаге. SSA-S выбор собственных троек и параметра Ь происходит один раз на втором шаге. В то время как в методе SSA-R выбор собственных троек требуется осуществлять на каждом шаге прогноза, что затрудняет использование данного метода в автоматическом режиме. Поэтому было принято решение в работе использовать метод прогноза SSA-S.

1.3.3 Прогнозирование ВР на основе метода группового учета аргументов

Основная идея прогнозирования ВР на основе метода группового учета аргумента (МГУА) (англ. - Group method of data handling, GMDH) базируется на гипотезе о том, что значение xN ВР Xi зависит от B его предыдущих значений:

XN = F(xN-1,XN-2, ... ,XN-в) , (1.13)

где в общем случае функция F оказывается неизвестной, однако, известна таблица ее значений:

XN XN-1

xn-1 xN-2 • • • xN-B XN-2 XN-3 • • • XN-B-1

(1.14)

хм—к хм—к—1 хм—к—2 • • • хм—В—к здесь N — В — к ^ 0.

Для решения подобных задач был разработан МГУА [36], согласно которому неизвестную функцию Г(хм—1 —2,... —В), связывающую входные переменные хм—1,хм—2, ■ ■ ■ ,хм—В и выходную переменную хм ищут в виде полинома Колмогорова-Габора [37]:

т т т

хм = ао + ^ ахм—г + ^ I] агзхм—¿хм—з +

¿=1 г=1 з=1 (1 15)

т т т

^2 ^2 аг,з,кхм—гхм—зхм—к,

г=1 з=1 к=1

где т ^ В.

Коэффициенты полинома Колмогорова-Габора в (1.15) вычисляются с использованием метода наименьших квадратов на основе данных, приведенных в (1.14). Из (1.15) видно, что точность аппроксимации функции Г (хм-1,хм-2, • • •, хм-В) в (1.13) зависит от порядка полинома Колмогорова-Габора и значения т. При увеличении числа т и порядка полинома Колмогорова-Габора сложность вычисления полинома Габора резко возрастает.

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

Также известен самоорганизующийся МГУА [101], который с вычислительной точки [38] зрения, оказывается менее затратным, нежели выше описанный классический МГУА. В данном методе для всех возможных пар сочетаний переменных {хм ,хм-1, •••, хм-В} - {х<п>,х<т>}, п = N - - 2, • • • N - В + 1, т = п - 1, ••• N - В строится множество моделей аппроксимирующих зависимость (1.15) первого уровня с помощью полиномов Колмогорова-Габора второго порядка:

х

(п,т) N

а,

(п,т)

+

а

(п,т)

х

(п)

+

а

(п,т)

х

(т)

+

а

(п,т) 1,1

х(п)х(т) •ЯУ гЛ/ .

(116)

где р = 1,2, ••• (В - 1)В/2.

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

р

р

р

р

р

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

1. Используя известные значения прогнозируемого ВР х0 = X(¿0),х1 = X(¿1),..., хм—1 = X—1), построить, используя классический или самоорганизующийся GMDH, GMDH-модель ВР (функцию Г в (1.13).

2. Последовательно вычислить прогнозируемые значения ВР:

Отметим, что в ряде случаев качество прогнозирования с помощью МГУА оказывается недостаточно высоким (например, при прогнозировании решений системы дифференциальных уравнений СДУ, описывающих аттрактор Лоренца [39]), поскольку здесь используется относительно небольшое число моделей-кандидатов, в связи с чем существует ненулевая вероятность потери значимой переменной в процессе отбора GMDH-моделей. Кроме того в ряде случаев могут возникать вычислительные проблемы в ходе решения систем уравнение метода наименьших квадратов [40; 41], а также известны примеры когда GMDH-модели ВР, порожденных ДС с относительно простой, с физической точки зрения, динамикой, оказываются весьма сложными [42; 43].

Один из активно используемых в настоящее время подходов к прогнозированию ВР основан на использовании искусственных нейронных сетей (ИНС) [27; 44; 45]. С точки зрения рассматриваемой задачи ИНС, как и рассмотренные выше формальные методы прогнозирования ВР, представляет некоторую функцию Г(хм—1,... ,Хм—1—В), зависящую от В переменных (В - база прогноза), которая описывает связь между прогнозируемым значением хм ВР хг и значениями ВР {хм — Ь . . . ,Хм—1—В}:

хм = Р(хм—ихм—2, ...,.х„-В) х>м+1 = Г (х'м хм—1, ... хм+1—В).

(1.17)

1.3.4 Искусственные нейронные сети

хм = Г(хм — Ь . . . ,Хм—1—ВX

(118)

Xi

Xo

X

N

— y

Рисунок 1.11 — Структура нейрона

Входной

Скрытый Выходной

слой

xn-1 xn-2

Xn—B

слой

Xn

Рисунок 1.12 — Структура простой ИНС

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

Отметим, что наилучшие результаты в рассматриваемой задаче, как показывает анализ опыта применения ИНС, достигаются при использовании рекуррентных нейронных сетей (RNN) с долгой краткосрочной памятью (Long short-term memory, LSTM) [46], с помощью которых, в ряде случаях, удается получать прогнозы долгосрочных зависимостей. Пример искусственного нейрона представлен на рисунке 1.11. Типичная архитектура ИНС представлена на рисунке 1.12.

Алгоритм прогноза с помощью ИНС реализуется выполнением следующей последовательности действий:

1. На основе результатов анализа прогнозируемого ВР x0 = X(t0)x = X(t1),..., xn-1 = X(tN-1), выбрать архитектуру ИНС и ее параметры:

- число нейронов Nneuron,

- функцию активации нейрона f,

- функцию инициализации начальных весов ИНС fw,

- функцию оптимизации ИНС fopt,

- число эпох обучения ИНС ^роась,

- число входов ИНС ^прЫ (база прогноза В).

- размер обучающей выборки г,

- размер тестовой выборки

- число выходов ИНС Хгпри (дальность прогноза Я).

2. Обучение ИНС.

3. Вычислить последовательно прогнозируемые значения ВР согласно 1.18.

Однако, оказывается, что ИНС с архитектурами NAR-net и LSTM-net,

используемые для прогнозирования солнечной активности, не обеспечивают требуемой точности долгосрочного прогноза [47; 48].

Недостатком применения ИНС в задаче прогнозирования ВР, порожденных ДС, описываемых с помощью феноменологических моделей, следует отнести отсутствие научно обоснованных рекомендаций, основанных на результатах анализа свойств прогнозируемого ВР, по выбору архитектуры ИНС, перечисленных выше параметров ИНС и алгоритмов ее обучения [49—52].

1.3.5 Метод коррекции прогнозов при использовании формальных методах

прогнозирования ВР

При использовании формальных методов прогнозирования ВР можно использовать следующий очевидный метод коррекции спрогнозированных значений, блок-схем0а которого представлена на рисунке 1.13, из которого видно, что алгоритм коррекции реализуется выполнением следующей последовательности действий:

1. Выбор метода прогнозирования ВР хг, где г = 0^ — 1 из множества < / >=< ЛЯ,ЗБЛ,СМВН,ЬБТМ >.

2. Идентификация параметров прогностической модели на основе использования значений г = N — В — 1, ..., N — 1-го значений прогнозируемого ВР -вектора базы прогноза хВ.

3. Прогнозирование значения N-го члена ВР - и его использование для предсказания дальнейшего поведения ДС.

4. Получение наблюдаемого значения хмоЪзег>.

5. Обновление базы прогноза хм_В_ 1 = хм—В,... ,хм—1 = хмоЪзег>.

6. Перейти к п. 2.

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

1.3.6 Итоги анализа формальных методов прогнозирования

1. Формальные методы прогнозирования ВР обеспечивают возможность вычисления рекуррентного прогноза ВР г = — 1 согласно (1.17) ВР, порожденных ДС, немеющими математических моделей.

2. Для коррекции рекуррентных прогнозов значений ВР х/м+1, х/м+2,..., вычисленных с помощью формальной модели ВР, порожденного ДС, описывающихся с помощью феноменологических моделей, на основе использования наблюдаемого значения ВР х°мзег необходима повторная идентификация значений параметров формальной модели ВР (при использовании ИНС - ее переобучение) на основе ВР {х2,х3,... ,хм—1,х°^} вычисление новых прогнозируемых значений ВР х^+1, х/м+2,... согласно (1.17).

3. Представляет практический интерес разработка новых подходов к коррекции прогнозированных значений ВР, порожденных ДС, описывающихся феноменологическими моделями, в которых не требуется повторная идентификация параметров прогностических моделей ВР при получении каждого следующего наблюдаемого значения ВР.

Рисунок 1.13 — Схема алгоритма коррекции прогнозов при использовании формальных моделей ДС, породивших прогнозируемый ВР

1.4 Классический метод коррекции прогноза Data Assimilation

Рассмотрим известный метод коррекции прогнозируемых значений ВР, порожденных ДС, описывающихся математическими моделями, называемый методом ассимиляции данных (англ. - Data Assimilation).

Определение 8. Data Assimilation (DA) - метод коррекции прогноза состояния ДС, данного на основе использования ее известной математической модели, с учетом данных наблюдений фактического состояния ДС.

Напомним, что методы коррекции расчетов, полученных на основе использования математических моделей ДС, данными наблюдений этих систем, получившие название методы ассимиляции данных, впервые были предложены для решения задач математической геофизики и предсказания погоды [53; 54]. Заметный всплеск интереса к данному методу, наблюдающийся в последние годы, связан, в первую очередь, со значительным прогрессом в развитии средств вычислительной техники и численных методов. Наибольшее распространение методы DA получили в метеорологии и физической океанографии, а также в оперативном численном предсказании атмосферных и океанских течений [55; 56]. В настоящее время теоретические и практические идеи ассимиляции данных применяются для решения различных инженерных (в том числе, разработке автоматизированных систем управления техническими системами) [57—59], математических [60—63] и геофизических задач [64—67].

В методе DA полагается, что эволюция исследуемой ДС описывается в дискретном времени следующей системой уравнений [54]:

Xk+i = M(xk,tk) + Wk, (i 19)

yk = H(Xk ,tk) + Vk,

где X - вектор состояния динамической системы, M - оператор или функция перехода, определяющая эволюцию системы во времени, Xk - вектор состояния системы в момент tk, xk+1 - вектор состояния системы в момент времени tk+1, yk - вектор наблюдений, H[Xk,tk) - оператор наблюдения, связывающий многомерное состояние системы с одномерным вектором наблюдений, wk, vk - ошибки модели и наблюдения, соответственно.

Алгоритм коррекции прогноза DA состоит из следующих этапов:

k k+1 k+2 k+3 t

Рисунок 1.14 — Схема работы метода Data Assimilation

1. Вычисление прогнозируемого значения вектора состояния ДС, первое уравнение в (1.19), называется уравнением прогноза.

2. Сравнение спрогнозированного вектора состояний ДС и наблюдаемого вектора состояния ДС, на основе применения второго уравнения в (1.19), называется уравнением наблюдения.

3. Корректировка вектора состояния ДС и вычисление следующего прогнозируемого значения ДС на основе использования скорректированного вектора состояния ДС.

Схема алгоритма DA представлена на рисунке 1.14.

Метод DA был разработан для линейной ДС, в математической модели в которой операторы M и H в (1.19) являются матрицами, поэтому:

M (xk ,tk ) = M xk( t) = M xk = Mk,

H (Xk ,tk ) = H Xk( tk) = H Xk = Hk,

а wk, vk - выборки из нормального белого шума с известными значениями ковариационных матриц ошибок Qk,Rk.

Из рисунка (1.14) видно, что вектор состояния ДС строится на основе спро-

_ f

гнозированного вектора состояния ДС xk, и далее на его основе вычисляется скорректированное значение XI. При этом учитывается значение наблюдения вектора состояния ДС yk только в текущий момент времени tk :

xa = Xfk + K y - H Xk )), (1.20)

где К - коэффициент передачи, е - расчетная ошибка (отклонение спрогнозированного на основе математической модели ДС вектора наблюдений ДС от реального значения вектора наблюдений). Метод DA призван минимизировать

к

ошибку е = у к — Н (х{), что в линейном случае выглядит как:

е = ук — Н Хк = Ук — Нк. Также на этом шаге рассчитывается скорректированная матрица ковариа-

а к,

ции Раа, согласно

Рак = (I — Кк Нк) Рк, (1.21)

здесь I - единичная матрица.

В предположении о том, что в (1.19) М и Н линейные операторы, математическое ожидание ошибки равняется нулю, а сами шумы и V к являются выборками нормального белого шума с ковариационными матрицами состояния Як и наблюдения Як, коэффициент передачи К рассчитывается по формуле:

Кк = Р{ НТ (Нк р{ НТ + Як)—1, (1.22)

£

где Рк - ковариационная матрица ошибок прогноза, вычисляется согласно формуле:

р{ = Мк (I — К—Н— )р к—МТ + як. (1.23)

На этапе корректировки прогноза с учетом вектора наблюдений вычисляют уточненное значение вектора состояния ДС по формуле:

ха = х{ + Кк (Ук — Н(х{)) = х{ + Кк (Ук — Нк), (1.24)

а также новую ковариационную матрицу ошибок прогноза:

Р к+1 = М+1 + Як, (1.25)

_ £

используемую далее для вычисления Хк+1 в будущий момент времени £к+1 согласно (1.19):

Х к+1 = М (ха). (1.26)

Затем повторяется описанная выше последовательность действий.

Предсказание эволюции ДС требует использование дополнительной информации, необходимой для выбора начальных условий и выбора параметров математической модели в (1.19), которую можно получить с помощью наблюдений. В этой связи возникает задача нахождения при заданной функции наблюдений yk (t) неизвестных параметров математической модели ДС (например, начального условия) таких, чтобы вектор состояния x удовлетворял (1.19), а вектор H был близок в некотором смысле к yk(t). Найденное решение x называется оценкой вектора состояния ДС и обозначается xa [54].

Далее по получению измеренных значений вектора состояния ДС рассчитывается скорректированный прогноз 1.20 и обновляется ковариационная матрица ошибок 1.21 Отметим, что коррекция прогноза в алгоритме DA с целью уменьшения вычислительных затрат может производиться не на каждом шаге по времени, но через несколько шагов по времени, что оказывается весьма важным для многомерных ДС.

Одним из популярных фильтров, лежащим в основе метода DA является фильтр Калмана, который использует статистическую интерполяцию. На основе стандартного фильтра Калмана разработано большое число модификаций, в настоящее время на практике используются различные модификации метода DA, в том числе:

- The Ensemble Kalman Filter (EnKF) - ансамблевый фильтр Калмана;

- The Stochastic Ensemble Kalman Filter (SEnKF) - стохастический ансамблевый фильтр Калмана;

- The Extended Kalman Filter (EKF) - расширенный фильтр Калмана нелинейный;

- The Unscented Kalman Filter (UKF) - нелинейный фильтр Калмана;

- The Original Ensemble Square Root Filter (EnSRF) - ансамблевый фильтр Калмана с квадратным корнем;

- The Local Ensemble Transform Kalman Filter (LETKF) - локальный ансамблевый фильтр Калмана;

- The Error-Subspace Transform Kalman Filter (ESTKF) - фильтр Калмана с преобразованием подпространства ошибок.

Далее рассмотрим основные фильтры, применяемые в методике DA, которые обеспечивают приемлемые показатели точности прогноза.

1.4.1 Фильтр Калмана

Одним из наиболее распространенных фильтров, лежащих в основе DA, является фильтр Калмана, названный в честь его автора Р. Калмана [54]. Фильтр Калмана (The Kalman Filter, KF) - это рекурсивный фильтр, с помощью которого оценивают вектор состояния ДС. При вычислении оценок вектора состояния ДС используется некоторый ряд наблюдаемых значений вектора состояния ДС, в которых присутствует шум.

Фильтр Калмана в отличие от других фильтров, схожих с ним по принципу работы, оперирует не только оценками вектора состояния, но и оценками его неопределенности плотности распределения вектора состояния, опираясь на формулу условной вероятности Байеса. Схема двухэтапного алгоритма DA на основе фильтра Калмана представлена на рисунке 1.15. (Здесь k - шага отсчета, m - длина прогноза.)

1. Этап прогноза (на этом этапе фильтра экстраполирует значения переменных состояния, а также их неопределенности).

2. Этап коррекции (на этом этапе фильтр по данным измерения, которые были получены с некоторой погрешностью, уточняет результат экстраполяции).

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

- во время моделирования ДС фоновые ошибки могут, как увеличиваться, так и уменьшаться;

- существует возможность выделения областей, в которых на предыдущем шаге анализа фоновые ошибки с помощью наблюдений удается уменьшить;

- существует возможность выделения областей, в которых из-за модельных ошибок фоновые ошибки выросли.

Начало

Цикл г от к до к + т

/

модели х0, исходном матрицы ковариации ошибок Р0

Вычисление матрицы Калмана К согласно (1.22)

Обновление оценки состояния х согласно (1.20)

Обновление матрицы ковариации ошибок по формуле (1.21)

Рассчитать скорректированный прогноз согласно (1.26)

Рассчитать новую ковариационную матрицу ошибок для следующего шага согласно Р/+1 = МшР%МГ+1 + Як

Увеличить г на 1

Вывести результат

Конец

Рисунок 1.15 — Схема алгоритма фильтра Калмана

1.4.2 Ансамблевый фильтр Калмана

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

Ансамблевый фильтр Калмана представляет собой вариант обобщённого фильтра Калмана, в котором ковариации ошибок прогноза оцениваются с помощью ансамбля прогнозов [69]. При этом применение такого подхода может быть реализовано с помощью распараллеливания вычислений с минимальными затратами на многопроцессорной вычислительной технике.

Различают два основных класса ансамблевых фильтров:

- стохастические;

- детерминированные.

Основное различие между этими схемами, основанными на ансамблевом подходе, состоит в способе генерации ансамбля исходных состояний для оценки [54].

Стохастические ансамблевые фильтры Калмана основаны на возмущении наблюдений - добавлении к истинным (реальным) значениям некоторого случайного шума с известной ковариационной матрицей.

При использовании детерминированных ансамблевых фильтров Калмана этап анализа, на котором вычисляются средние значения ансамбля и ковариационной матрицы ошибок наблюдений выполняется один раз [54]. В каждый из вариантов Еп^ используются:

- ансамблевый прогноз;

- обновление или анализ ансамбля.

Шаг ансамблевого прогноза одинаков для всех типов ансамблевых фильтров Калмана. Здесь изменение во времени каждого из членов ансамбля реализуется согласно математической модели ДС.

Схема алгоритма Еп^ представлена на рисунке 1.16, где к - шага отсчета, т - длина прогноза. Из рисунка 1.16 видно, что обсуждаемый алгоритм реализуется выполнением следующей последовательности действий:

1. Инициализация: формируется N наблюдения возмущенной модели для каждого члена ансамбля. Начальный ансамбль Еп^ формируется из исторических знаний:

2. Прогнозирование: вычисление спрогнозированных значений для каждого члена ансамбля.

3. Корректировка спрогнозированных значений: вычисление коэффициента передачи К.

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

Ут^ = Н (Х^) + Ут,М.

(1.27)

ХМ = Хм + КЛт,Ы • 5. Повторение п. 1-3 для каждого шага прогноза.

(1.28)

Рисунок 1.16 — Схема алгоритма ансамблевого фильтра Калмана

Напомним, что существуют различные схемы ансамблевого фильтра Калмана: ансамблевый трансформированный фильтр Калмана ETKF, ансамблевый корректирующий фильтр EAKF, ансамблевый фильтр квадратного корня EnSRF и др. [53; 54; 70; 71]. Одним из усовершенствований фильтра EnKF является стохастический подход The Stochastic Ensemble Kalman Filter (SEnKF) в отличие от детерминированного подхода в ансамблевых фильтрах Калмана генерация

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

1.5 Постановка задач исследования

Анализ состояния предметной области, проведенный в данной главе, позволяет сделать следующие выводы:

1. Задача прогнозирования ВР, порожденных ДС, описываемых, как с помощью математической модели, так и с помощью феноменологической модели, является актуальной в различных областях науки и человеческой деятельности.

2. Методы прогнозирования ВР в зависимости от типа ДС, породившей данный ВР, можно разделить на два независимых класса:

- методы прогнозирования ВР, основанные на использовании математических моделей ДС, породившей изучаемый ВР;

- формальные методы прогнозирования ВР, основанные на использовании феноменологической модели ВР.

3. На практике известны случаи, когда требуется прогнозировать ВР, у которых оказываются известными как ранее спрогнозированные, так и соответствующие наблюдаемые значения, порождённые обоими типами ДС, что потенциально, обеспечивает возможность уточнения прогнозов значений ВР.

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

5. Для коррекции значений ВР, спрогнозированных на основе использования математической модели ДС, породившей данный ВР, и соответствующего наблюдаемого значения ВР, используют метод DA.

6. Не существует методики, обеспечивающей совместное использование формальных методов прогнозирования ВР и метода DA.

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

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

Задачи исследования:

1. Обоснование возможности интеграции формальных методов прогнозирования ВР и метод DA.

2. Разработка методики интеграции методов прогнозирования ВР и метода

DA.

3. Разработка программного обеспечения, поддерживающего методику интеграции формальных методов прогнозирования ВР и метода DA, а также известные методы прогнозирования ВР.

4. Сравнительный анализ качества прогнозирования ВР на основе использования формальных методов прогнозирования ВР и методики интеграции формальных методов прогнозирования ВР и метода DA.

Глава 2. Разработка методики интеграции формальных методов прогнозирования ВР и метода DA

2.1 Обоснование возможности интеграции формальных методов прогнозирования ВР и метода DA

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

Напомним, что математическая модель изучаемой ДС представляет собой дифференциальное уравнение второго порядка:

x + ш^х = 0, (2.1)

где шо = 2п/о - циклическая частота колебаний, /о - частота колебаний, с начальными условиями х(0) = х0, Х(0) = v(0), а его решение в наиболее общем виде есть функция вида:

x(t) = A cos (w0t) + B sin (w0t), (2.2)