Abstract

Abstract

The paper proposes a method of short-term forecasting based on the identification of the spanning algebraic sequence c using evolutionary algorithms. This method is especially effective in the case of short time series. When there is not enough data for training standard models, the developed method still allows to extract the maximum available information on the process behavior in order to extrapolate it to future points in time.

Full Text

Введение Прогнозирование временных рядов является одной из наиболее актуальных проблем, возникающих в науке и технике. Условно говоря, временных рядов можно разделить на долгосрочное прогнозирование и краткосрочное прогнозирование. Целью долгосрочного прогнозирования временных рядов является построение модели процесса, которая затем используется для экстраполяции поведения в прошлом на достаточно отдаленное будущее. Цель краткосрочного прогнозирования временных рядов, в общем-то, остается той же, т.е. также строится модель и предсказывается поведение в будущем. К сожалению, методы, разработанные для моделей долговременных временных рядов, часто непригодны для кратковременных временных рядов по причине недостаточного объема имеющихся данных. С другой стороны, при краткосрочном прогнозировании временных рядов достаточно лишь иметь возможность делать прогноз на один шаг вперед (иногда сопровождая это поиском локальных минимумов и максимумов). Методы долгосрочного и краткосрочного прогнозирования временных рядов основаны на различных подходах и моделях. Среди них можно выделить различные техники сглаживания, такие как фильтр скользящего среднего и экспоненциальное сглаживание [1-4], методы, основанные на использовании модели типа авторегрессии проинтегрированного скользящего среднего (АРПСС) [5-8], методы искусственного интеллекта типа НС (искусственная нейронная сеть) [9-12] и др. 1. Стратегия прогнозирования Предположим, что для построения прогнозной модели процесса имеется наблюдение: (1) где - наблюдение для текущего момента времени. На основе нечеткого числа наблюдений можно построить минор Ганкеля . Пусть (если тогда последовательность будет являться алгебраической прогрессией, поэтому нахождение значения не составляет труда). Сделаем еще одно предположение. Пусть последовательность (1) является результатом сложения алгебраической последовательности и случайного шума. Другими словами, считается, что где - аддитивный шум, и , где . Более того, будем также предполагать, что бесконечная последовательность является алгебраической прогрессией, которая выступает в роли остовной последовательности, определяющей глобальную динамику временного ряда. Задача состоит в том, чтобы каким-то образом выделить аддитивный шум Используя простые рассуждения, можно показать, что у данной задачи существует бесконечное число решений. Поэтому будем искать такой аддитивный шум, который бы вносил минимальные искажения в исходный временной ряд. Таким образом, основной задачей будет являться нахождение поправок которые максимизируют следующую целевую функцию: (2) где , Очевидно, что ; . Если целевая функция достигает максимума при причем Параметр a устанавливает относительные значения штрафов, накладываемых на целевую функцию величиной определителя и взвешенной суммой поправок (в случае оба значения штрафа будут иметь одинаковый вес). Коэффициенты формирует коридор допустимых значений поправок Все поправки будут иметь одинаковый вес, если Чем больше b, тем выше вес поправки для текущего наблюдения по сравнению с более ранними наблюдениями. Другими словами, коридор допустимых значений поправки в текущий момент времени ýже, чем в более ранние моменты времени. Такое поведение желательно, если используется предположение о том, что значимость наблюдения зависит от его близости к текущему моменту времени. Для нахождения почти оптимального набора поправок использовалcя эволюционный алгоритм (ЭА). Настройка параметров алгоритма осуществлялась путем тестирования на примере синтетического временного ряда. Для начала генерировалась периодическая последовательность (значения ее семи элементов за один период равнялись 0,5; 0,7; 0,1; 0,9; 0,3; 0,2; 0,8). Эта последовательность образует остовную алгебраическую последовательность. Далее ко всем элементам последовательности добавлялись случайные числа, равномерно распределенные в интервале (рис. 1). Полученный таким образом тестовый временной ряд использовался для проверки работоспособности предложенного метода прогнозирования. Рис. 1. Тестовый временной ряд (сплошная линия) и образующая его периодическая последовательность (пунктирная линия) Рис. 2. Связь между абсолютным значением определителя матрицы Ганкеля и ее порядком Первым шагом будет являться выделение базового фрагмента временного ряда, который требуется для восстановления остовной алгебраической последовательности. Для заданной периодической последовательности решение этой задачи не составляет труда (ее H-ранг будет равен 8), однако для временного ряда задача усложняется из-за наличия аддитивного шума в данных. На рис. 2 показана зависимость абсолютного значения определителя матрицы Ганкеля от размерности матрицы. Видно, что определитель не равен нулю даже в случае, когда размерность матрицы равна 8. Тем не менее мы зафиксировали размерность матрицы Ганкеля на уровне 8 длина базового фрагмента временного ряда Оценка качества прогнозирования временного ряда на рис. 1 для матриц Ганкеля других размерностей приведена в последующих параграфах. 2. Использование эволюционных алгоритмов Естественным образом возникает вопрос, можно ли использовать иной метод оптимизации для поиска почти оптимального набора поправок . Предложено использовать эволюционные алгоритмы (ЭА). Хотя применение ЭА не всегда гарантирует наилучшее решение всех проблем, неоспоримыми преимуществами данных методов являются эффективная работа с большим количеством переменных, большим объемом экспериментальных данных и нелинейными аналитическими функциями, а также пригодность для оптимизации целевых функций с особенно сложной топологией [13]. Хромосомы особей, которые используются для кодирования набора поправок, включают генов, принимающих вещественные значения. Начальная популяция состоит из m особей (в нашем случае значения генов которых генерируются случайным образом в интервале . Функция приспособленности каждой особи определяется уравнением (2). Для улучшения начальной популяции производится отбор (селекция) четного количества особей. В качестве метода селекции использовался случайный метод рулетки [14]. Чем выше уровень приспособленности особи, тем больше вероятность, что ее выберут для последующего скрещивания. Однако вероятность выбора особи с малым значением функции приспособленности не равна нулю. Кроме того, в результате селекции могут быть выбраны одинаковые копии (клоны) одной и той же особи. Далее все особи разбиваются по парам. Скрещивание особей производится для всех образованных пар. Для этого использовался одноточечный модифицированный оператор скрещивания b-кроссовер [15] (положение точки случайно для каждой из пар). При этом происходит не только обмен генами, но также стремятся к тому, чтобы два новых потомка были более или менее схожи друг с другом, а также с родительскими особями. Алгоритм b-кроссовера описывается следующими уравнениями: где k - номер поколения; - j-й ген хромосомы левого родителя в паре; - j-й ген хромосомы правого родителя в паре; - j-й ген хромосомы левого потомка; - j-й ген хромосомы правого потомка; минус в нижнем индексе означает, что j-й ген располагается выше точки скрещивания; плюс в нижнем индексе означает, что j-й ген располагается ниже точки скрещивания. Во всех случаях j-й ген хромосомы потомка лежит в интервале . Оператор b-кроссовера сближается с классическим алгоритмом скрещивания, когда b стремится к нулю или бесконечности (если , то оба потомка будут иметь одинаковые хромосомы). Обычно оператор скрещивания применяется не ко всем парам особей в промежуточной популяции [16]. Вероятность того, что пара особей будет участвовать в скрещивании, определяется коэффициентом k. Для предотвращения сходимости к локальному максимуму и реализации глобальной стратегии поиска решений используется оператор мутаций. Интенсивность мутационного процесса определяется параметром m [13]. Алгоритм заключается в замене генов всех особей в текущей популяции на случайную величину, равномерно распределенную в интервале . 3. Выбор параметров эволюционного алгоритма В общем случае выбор параметров ЭА осуществляется эмпирически, хотя в [17, 18] описываются некоторые известные принципы их выбора. Перед применением ЭА исследователь должен заранее установить значения для следующих параметров: коэффициента скрещивания k, параметра мутаций m, коэффициента сходства b и числа поколений. Мы будем следовать рекомендациям для классической модели ЭА [14]. Коэффициент скрещивания k будет выбираться из интервала , параметр мутаций m - из интервала . Выше уже отмечалось, что в классической модели ЭА коэффициент сходства (или Однако мы предпочитаем, чтобы значения данного коэффициента находились в интервале В литературе нет четких указаний относительно требуемого числа поколений; нередко для хорошей сходимости требуется, чтобы их было не менее 80. В данной работе при проведении вычислительных экспериментов число поколений равнялось 40. За одно выполнение ЭА формируется один набор поправок . Очевидно, что результат зависит от начальной популяции особей (помимо других случайных факторов). Поэтому выполнение ЭА (при фиксированных значениях параметров) осуществлялось 100 раз, далее вычислялось 100 прогнозных оценок а затем определялась среднеквадратическая ошибка прогноза. Для достижения лучшей точности прогноза значения параметров k, m и b варьировались (табл. 2). Наилучший результат получен при и В дальнейшем использовался данный набор значений параметров. Можно повысить эффективность метода прогнозирования на базе ЭА путем выбора соответствующего значения параметра в (2), который отвечает за коридор допустимых значений поправок, а также путем исключения выбросов в прогнозных оценках. Вычислительные эксперименты показали, что наилучшие результаты получаются при (таблица, жирный шрифт) и исключении 10 % выбросов (рис. 3). Рис. 3. Прогнозирование тестового временного ряда с использованием ЭА для нахождения поправок : а - без удаления выбросов; б - удаляется 10 % выбросов от общего числа оценок; в - удаляется 25 % выбросов от общего числа оценок Среднеквадратические ошибки прогноза на примере тестового временного ряда Метод b 0 0,25 0,5 0,75 1 1,25 1,5 1,75 2 3 4 ЭА, β = 0 0,1899 0,1916 0,1919 0,1939 0,1878 0,1883 0,1886 0,1921 0,1948 0,1904 0,2061 ЭА, β = 2 0,1870 0,1885 0,1886 0,1883 0,1861 0,1886 0,1893 0,1882 0,1896 0,1888 0,1893 4. Вычислительные эксперименты на примере реального временного ряда Для апробации предложенного подхода на примере реального временного ряда использовались суточные замеры температуры по сухому термометру, приведенные в [19]. Первоначальная задача состоит в идентификации базового фрагмента временного ряда. Первый минимум определителя матрицы Ганкеля достигается, когда ее размерность равна 11 ( соответственно, длина базового фрагмента временного ряда ). Полученная длина базового фрагмента фиксировалась для всего временного ряда. Значения параметров ЭА были такие же, как и в случае тестового временного ряда. Прогноз осуществлялся последовательно на один шаг вперед (рис. 4). Нетрудно заметить, что исключение 10 % выбросов (как и для тестового временного ряда) улучшает точность прогнозирования. Для сравнения с другими методами прогнозирования использовался тот же самый временной ряд. Прогноз осуществлялся по методу анализа временных рядов Бокса - Дженкинса на базе модели АРПСС (3,0,3) [20], так как экспериментально было установлено, что такая структура модели, 3-0-3 (т.е. когда порядок авторегрессионной составляющей равен 3, порядок интегральной составляющей равен 0, а порядок составляющей скользящего среднего равен 3), лучше всего описывает имеющиеся данные. Среднеквадратическая ошибка прогноза составила 0,455 6, что почти в два раза меньше значения, полученного по нашему методу. На рис. 5 приведены результаты прогноза временного ряда для обоих методов (для наглядности прогнозирование на базе АРПСС начиналось с 21-го элемента временного ряда). Рис. 4. Прогнозирование реального временного ряда с использованием ЭА для нахождения поправок : а - без удаления выбросов; б - удаляется 10 % выбросов от общего числа оценок; в - удаляется 25 % выбросов от общего числа оценок На первый взгляд, модель АРПСС (3, 0, 3) имеет преимущество перед предложенным методом. Однако на рис. 5 видно, что модель АРПСС демонстрирует поведение, характерное для процесса скользящего среднего. Фактически все острые пики исходной кривой чрезмерно сглаживаются из-за происходящего в модели АРПСС (3, 0, 3) усреднения входных данных. С другой стороны, предложенный метод существенно отличается от известных статистических алгоритмов, поскольку основывается на локальной идентификации остовной алгебраической последовательности на каждом временном шаге независимо от результатов, полученных на предыдущих шагах. Тем не менее следует признать, что в некоторых точках прогнозные оценки по предложенному методу значительно отличаются от истинных значений временного ряда, что в конечном счете ухудшает качество прогноза среднего значения. В то же время предложенный метод хорошо распознает положение локальных флуктуаций исходного временного ряда. Данное свойство позволит, например, лучше предсказывать локальные максимумы и минимумы на сутки вперед, чем процесс АРПСС [2]. Рис. 5. Прогнозирование реального временного ряда: а - на базе предложенного метода с использованием ЭА для нахождения поправок; б - на базе АРПСС (3, 0, 3); в - на базе смешанного метода Интересным развитием предложенного подхода является объединение процесса АРПСС (3, 0, 3), который обладает свойствами скользящего среднего, и нашего метода, который лучше воспроизводит локальные изменчивости временного ряда. На рис. 5 представлены результаты смешанного метода, в котором вычисляется среднее арифметическое значение прогнозных оценок, полученных обоими методами. Хотя среднеквадратическая ошибка такого метода выше, чем у АРПСС (3, 0, 3), однако вид получаемой кривой лучше соответствует исходному временному ряду. Очевидно, что у любого метода есть своя специальная область применения, в которой он проявляет себя наилучшим образом. Предложенный метод базируется на предположении о том, что в основе каждого процесса, описываемого временным рядом, лежит некоторый детерминированный закон его развития во времени. Ясно, что метод будет работать уже не так хорошо, если элемент случайности для временного ряда превалирует над детерминированной моделью процесса. Кроме того, предложенный метод существенно снижает свою эффективность, когда длительность имеющегося временного ряда меньше, чем минимально допустимая длина остовной алгебраической последовательности, определяющей модель процесса. Возможным направлением повышения точности прогноза в предложенном методе является дополнительное совершенствование алгоритма идентификации остовной последовательности. Более того, на данный момент длина базового фрагмента временного ряда является фиксированной величиной. Поэтому еще одно допущение, лежащее в основе предложенного метода, касается стационарности процесса. Возможность изменения длины базового фрагмента в зависимости положения точки начала прогноза временного ряда, а также совершенствование алгоритма идентификации остовной последовательности является темой дальнейших исследований. Заключение Разработана и протестирована новая методика краткосрочного прогнозирования временных рядов. Она основана на выявлении остовной алгебраической последовательности в имеющихся данных наблюдений. Принимая во внимание, что далеко не всегда достаточно данных для идентификации математических моделей, описывающих эволюцию временных рядов, предложен подход определения остовной алгебраической последовательности из короткого временного ряда, позволяющий извлечь как можно больше информации о модели процесса, насколько это возможно, а затем использовать эту модель для экстраполяции прошлого поведения в будущем. Предложенный метод прогнозирования базируется на предположении о том, что в основе каждого процесса, описываемого временным рядом, лежит некоторый детерминированный закон его развития во времени. Таким образом, представленный метод прогнозирования основан на еще одном предположении: что процесс стационарен. Тем не менее это важный прогресс в области прогнозирования, поскольку не существует иных методов прогнозирования, которые могут справиться с краткосрочными случайными процессами. Для идентификации остовной алгебраической последовательности эффективно применены эволюционные алгоритмы. Оригинальные функции пригодности вводятся для генетических алгоритмов с целью гарантировать оптимальную сходимость к почти оптимальным решениям. В итоге целью всех методов прогнозирования временных рядов является построение модели процесса, а затем использование этой модели для последних значений временного ряда для экстраполяции прошлого поведения в будущем. Такой подход хорошо работает для длинных временных рядов, где количество имеющихся данных позволяет строить такие методы прогнозирования, которые используют сотни или даже тысячи точек данных для построения математической (обычно аппроксимирующей) модели процесса. Но ситуация совершенно иная, когда доступное количество данных мало - тогда лучшим решением задачи является предлагаемый метод прогнозирования.