CC BY
40
убывания 1/k (а коэффициенты ak вообще не убывают, они и так малы - имеют величину порядка e), т.е. эта часть ряда сходится медленно и поэтому в решении необходимо учитывать все ее члены.
Подобными способами можно конструировать непрерывные замещающие функции для периодического сигнала (воздействия), имеющего в одном периоде несколько точек разрыва первого рода и определять необходимое число n членов ряда Фурье для обеспечения заданной точности вычислений при суммировании. При этом надо просуммировать все члены ряда, для которых ke остается малым, и в том случае, когда их число превышает n, полученное из оценки (1). Таким образом, чем меньше e и, следовательно, чем меньше замещающая функция отличается от исходной, тем большее число членов ряда необходимо суммировать.
Библиографический список
1. Смирнов В.И. Курс высшей математики: учеб. для вузов. М.: Наука, 1974. Т. 2. 656 с.
2. Фихтенгольц Г. М. Курс дифференциального и интегрального исчисления. М.: Наука, 1970. Т. 3. 656 с.
УДК 004.94
МЕТОД СТОХАСТИЧЕСКОГО КВАЗИГРАДИЕНТА В ЗАДАЧЕ АДАПТАЦИИ ПРОГНОЗИРУЮЩИХ МОДЕЛЕЙ
© И.А. Серышева1, Ю.П. Хрусталёв2
Иркутский государственный технический университет, 664074, Россия, г. Иркутск, ул. Лермонтова, 83.
Рассматриваются проблемы адаптации прогнозирующих моделей временных рядов (моделей АРСС), которые используются в задачах оценивания состояния динамических систем. Необходимость адаптации вызвана, в частности, ограниченной длиной временных рядов, по которым находятся первоначальные оценки параметров моделей. Предлагается для построения адаптивных процедур использовать метод обобщенного стохастического градиента. Показано, что применение данного метода обеспечивает хорошую сходимость процессов подстройки параметров.
Ил. 5. Библиогр. 10 назв.
Ключевые слова: недоопределенные системы; идентификация; прогнозирующие модели; процессы авторегрессии; критерии адекватности; временные ряды; оптимизация.
METHOD OF STOCHASTIC QUASIGRADIENT IN THE PROBLEM OF PREDICTIVE MODEL ADAPTATION I.A. Serysheva, Yu.P. Khrustalev
Irkutsk State Technical University, 83 Lermontov St., Irkutsk, 664074, Russia.
The paper considers the adaptation problems of time series predictive models (autoregression moving average models) that are used in the estimation problems of dynamical system state. The need for adaptation is caused, in particular, by the limited length of the time series, which are used for finding the initial estimates of model parameters. It is proposed to use the method of generalized stochastic gradient for building adaptive procedures. The application of this method is shown to ensure good convergence of parameter adjustment processes. 5 figures. 10 sources.
Key words: underdetermined systems; identification; predictive models; autoregression processes; adequacy criteria; time series; optimization.
Практически все измерительные процедуры основаны на операции сравнения измеряемой величины с некоторой эталонной величиной или мерой. При этом в соответствии с основным положением метрологии, всегда имеет место погрешность измерения. В результате мы можем в процессе обработки измерительной информации найти некоторую оценку измеряемой величины. Как правило, речь идет о нахождении в некотором смысле оптимальных оценок.
Для обеспечения единства измерений физических величин в стране создана система эталонов и существуют так называемые поверочные схемы, определяющие способы доведения точностных характеристик эталонов до потребителей. Верхним звеном поверочной схемы являются эталоны (первичный и вторичные). Таким образом, точность измерения той или иной физической величины определяется точностью эталонов.
1Серышева Ирина Анатольевна, старший преподаватель кафедры автоматизированных систем, тел.: (3952) 405164, e-mail: [email protected]
Serysheva Irina, Senior Lecturer of the Department of Automated Systems, tel.: (3952) 405164, e-mail: [email protected]
2Хрусталёв Юрий Петрович, кандидат технических наук, доцент кафедры вычислительной техники, тел.: (3952) 405107, email: [email protected]
Khrustalev Yuri, Candidate of technical sciences, Associate Professor of the Department of Computing Machinery, tel.: (3952) 405107, e-mail: [email protected]
Для повышения точности и надежности воспроизведения физических величин создаются групповые эталоны, элементами которых являются меры. В процессе эксплуатации эталонов величины, воспроизводимые мерами периодически сравниваются (сличаются) между собой. Задача обработки данных, получаемых в процессе функционирования групповых эталонов, заключается в нахождении оценок величин, воспроизводимых мерами по измеренным разностям этих величин.
Можно рассматривать различные схемы измерений:
- схему сличения всех элементов, входящих в эталон, с каким-то одним из них, выбранным в качестве опорного;
- схему сличения «каждый с каждым»;
- схему, основанную на разбиении всех элементов на отдельные подгруппы и измерениями, выполненными сначала внутри подгрупп, а затем - между подгруппами и т.д.
С точки зрения математики все эти схемы эквивалентны. Поясним выше сказанное. Пусть в групповой эталон входит п элементов. Тогда, при использовании схемы с опорным элементом имеется п - 1 результат измерений.
На рис. 1 представлена упрощенная схема системы измерений в групповом эталоне.
Рис. 1. Упрощенная схема измерений в групповом эталоне
Введем в рассмотрение следующие векторные величины:
Y - вектор состояния группового эталона; yfe - fc-ая составляющая вектора Y, fc = 1 ,2 ,. . ,,n; = у0 п - у;- вектор измерений, где у0п - значение физической величины, воспроизводящейся опорным элементом. Пусть это будет элемент т.
Тогда операцию измерений можно записать как
Z = Л ■ Y, (1)
где А - матрица наблюдений (будем использовать термин «наблюдение» как измерение, выполненное в отсутствие погрешностей [1]. Более точно как измерение, при котором уровень шумов существенно ниже уровня полезного сигнала и им можно пренебречь).
Матрица А имеет единицы в столбце с номером т, и «-1» в элементах матрицы А, соответствующих составляющей группового эталона. Для простоты будем считать опорным первый элемент, тогда первый столбец матрицы А содержит единицы. Диагональ, расположенная выше главной, состоит из -1. Все
остальные элементы матрицы наблюдений равны нулю.
Размерность матрицы наблюдений равна (п -1 ) х п. Таким образом, система (1) - недоопределен-ная и имеет бесчисленное множество решений. Единственное решение можно получить, вводя в систему какое-либо дополнительное ограничение. Простейшим ограничением является требование поиска решения с минимальной нормой вектора У (У - оценка вектора У). Тогда оценка опорного элемента уоп может быть найдена с помощью псевдослучайной матрицы £ [2].
Очевидно, что при отсутствии шумов измерительной системы достаточно найти оценку у; для какого-либо элемента группового эталона. Оценки всех остальных составляющих вектора находятся непосредственно из уравнения (1). Естественно выбрать в качестве такого элемента опорный элемент группового эталона.
Тогда справедливо [3]
Уоп = У = ^Г^1 г;. (2)
Заметим, что при использовании других схем измерений ранг матрицы не меняется, т.е. система остается недоопределенной. Усложнение измерительной схемы может быть полезной только при уменьшении уровня шумов, что, при принятом предположении об их малости, не дает никаких преимуществ.
Типичными представителями групповых эталонов являются эталоны времени и частоты, в которых измерения ведутся постоянно через определенные интервалы времени. Все выше сказанное остается справедливым для эталонов времени, с учетом того, что векторы У и 1 следует рассматривать как дискретные функции времени £, т.е. У = У( £), 1 = 1( £). Тогда эталон времени и частоты можно трактовать как динамическую систему, а процедуру нахождения оценок У( £) как задачу оценивания состояния динамической системы. В этом случае уместно рассматривать алгоритмы обработки измерительной информации, получаемой в процессе «ведения» эталонов как алгоритмы оценивания состояния динамических систем по результатам наблюдений [4]. Задача обработки данных в темпе поступления измерительной информации, т.е. в режиме реального времени, формулируется как задача фильтрации, решение которой в виде рекуррентных соотношений известно как фильтр Калмана. Поскольку в эталонах времени и частоты составляющие эталона (водородные стандарты) эксплуатируются в условиях, когда их состояния у - относительные отклонения частоты не зависят друг от друга, а связаны в единое целое в основном на уровне измерительной системы, целесообразно произвести декомпозицию фильтра. Тогда оценка опорного элемента в момент £ может быть получена из выражения [5].
уо п = &Е"= 1(г; + У;)' (3)
где - вес -го измерения, у; - прогноз I-ой составляющей.
В (3) для простоты не указан индекс £. Если считать все измерения равноточными, то & = Ц для всех I, формула (3) примет вид:
у0П=^Г= + (4)
В (3) и (4) для удобства использования «индексов» введено фиктивное измерение у0П - г0П = 0.
Анализ погрешностей алгоритмов, реализующих формулы вычисления значений (оценок) опорного элемента (2) и (4) показывает, что в последнем случае (т.е. в случае использования прогнозов) можно получить более точный результат.
Прогнозы у( £) можно вычислить на основе математических моделей, описывающих динамику системы. Такие модели являются динамическими стохастическими - моделями авторегрессии скользящего среднего (АРПСС). Методика построения моделей АРПСС хорошо известна. Имеются программные комплексы для анализа временных рядов, реализованные в пакетах прикладных программ STATISTICA и STAT-GRAPHICS.
Алгоритмы, реализующие методы, основанные на использовании прогнозирующих моделей (моделей АРПСС) при вычислении частоты (относительных отклонений частот) водородных стандартов, детально исследованы методами статистического моделирования и показали хорошие результаты при решении практических задач [6].
Однако при использовании таких алгоритмов требуется определенный опыт в построении моделей АРПСС по эмпирическим данным и соответствующая теоретическая подготовка специалистов. В реальных ситуациях эти требования выполняются далеко не всегда. Кроме того, для построения моделей используются временные ряды, длина которых ограничена. Во-первых, в рассмотрении исследователей, занимающихся построением моделей, имеются ограничения выборки, а во-вторых, при чрезмерном увеличении объема выборки можно столкнуться с нестационарностью процессов, что приведет к обесцениванию проделанной работы.
Естественным выходом в этой ситуации представляется использование адаптивных моделей временных рядов. В общем случае параметры моделей АРПСС представляются в виде вектора // = [Ф • 0 ] ,
где Ф = [0ф2.....фр ]т, 0 = [ 0, 02.....0, ]т, р - порядок авторегрессии (АР), - порядок скользящего среднего (СС). Как правило, для реальных временных рядов р и ц не превышают соответственно 3 и 2 [8]. В процессе адаптации (в отличие от самоорганизации) структура моделей (т.е. порядки р и ц) не меняется.
Цель адаптации - минимизировать сумму квадратов отклонений прогнозов частоты от их «истинных» значений. Прогнозы вычисляются как векторное произведение вектора параметров модели // = [Ф • 0] т на «вектор состояния процесса» ? = [у^ ^_ . ., ^_р •
] , где [ ] - вектор
коэффициентов АР, [ ] - вектор ко-
эффициентов СС, - ошибки прогноза. Таким обра-
зом, прогноз на один шаг вперед у( 1) равен
у( 1 ) = //т?=?т//. (5)
Целевая функция ( ) запишется при этом как
Я//)=£Г= 1(у(£)-у( 1 ))2, (6)
где - длина временного ряда (объем выборки).
Уидроу [10], рассматривая работу фазированной антенной решетки в режиме подстройки параметров, описывает ее с помощью дискретного адаптивного сумматора. При этом выходной сигнал он представляет как скалярное произведение вектора весовых коэффициентов на вектор входного сигнала составляющие которого задерживаются на тактов, т.е.
7 = (7)
Формула (7) полностью совпадает с выражением (5). Следовательно, результаты, полученные Б. Уидроу и его сотрудниками, можно использовать при решении нашей проблемы. В частности, авторы рассматриваемой работы называют среднее значение квадрата ошибки характеристической функцией адаптивного фильтра, а проблему адаптации весовых коэффициентов в стационарном режиме формулируют как оптимизационную задачу. В нестационарном случае процесс адаптации, авторы цитируемой работы трактуют, как слежение за низшей точкой впадины характеристической функции, положение которой может меняться во времени.
Рассмотрим сначала стационарный режим работы адаптивного фильтра. В нашем случае, т.е. при использовании моделей АРСС, такую ситуацию можно толковать как процесс подстройки параметров прогнозирующей модели, коэффициенты ф; и 0,, которой оценивались по ограниченной выборке. То есть речь может идти о применении стохастических квазиградиентных методов [8].
В соответствии с методом, предложенным Ю.М. Ермольевым [8], будем рассматривать последовательность точек
//5+ 1 =7Г/К//5-р5Х5С), 5 = 0 ,1,. . . (8)
где - оператор проецирования вектора параметров // в область допустимых значений; 5 - номер шага; р5 - величина шага в заданном направлении; - нормирующий множитель; ^- случайный вектор, условное математическое схождение которого
М((7у. ,,у5) = а/у(/?5) + Ь5, 5 = 0 ,1 ,. . . (9) где - случайная величина; - некоторый вектор, зависящий от последовательности ( );
^у(//5) - обобщенный градиент целевой функции (6).
Если в (6) положить , то
ЖС/у0,.■ ■ ,у5) будет равно обобщенному стохастическому градиенту, т.е.
ЖГ/У°.....у5) = ^(//5).
Поскольку в нашем случае целевая функция ( ) гладкая, то £у(/5) - стохастический градиент, найдется как производная от выражения (6), т.е.
^(//5) = [Я= 1(у(£)-У( 1 ))2 ]'. (10)
Учитывая (5) и обозначив ошибку прогноза как
= а(£)= у (£)-&( 1 ), (11)
и рассматривая процедуру адаптации только на одном
шаге (т.е. N = 1 и, стало быть, в (10) можно опустить знак суммы), получим алгоритм адаптации вектора параметров на одном шаге:
/+1 = 7^-0^-2 а5-/?5). (12)
Обозначив удвоенный множитель как у, т.е. у = 2 у5 и положив его равным 1 окончательно получаем
/+^(Г-рл-/) (13)
Таким образом, коррекция вектора параметров // прогнозирующей модели на очередном шаге пропорциональна ошибке прогноза.
Коэффициент р выбирается исходя из сходимости процесса адаптации. В работе Ю.М. Ермольева [8] показано, что процесс адаптации сходится с вероятностью 1 при выполнении условия
£Г=оР* = °°. (14)
Выполнение условия (14) обеспечивает сходимость процесса адаптации при любом выборе начального вектора //о. Поскольку операция проектирования гарантирует попадание вектора параметров в пределы области устойчивости коэффициентов авторегрессии и в пределы области обратимости [7] коэффициентов скользящего среднего, то можно это условие ослабить, предполагая, что при 5 — со р5будет стремиться к нулю. Это крайне важно сделать, учитывая так называемый «шум градиента» [10]. Шум градиента имеет место, когда находясь в точке оптимума /?*, мы продолжаем процесс адаптации, «дергая» без надобности коэффициенты авторегрессии и скользящего среднего. В результате качество функционирования системы ухудшается. Далее в стационарном случае при неограниченном увеличении длины временного ряда, влияние текущей ошибки прогноза на поправку коэффициентов и должно в пределе стремиться к нулю.
Естественен выбор коэффициента адаптации как 1 //с. То есть мы полагаем, что влияние последующих членов выборки (временного ряда) убывает в обратно пропорциональной зависимости от . Тогда на первом шаге адаптации этот коэффициент равен 1 (р 1 = 1).
И, наконец, если мы хотим обезопасить себя от небольших возмущений, действующих на параметры моделей, необходимо зафиксировать коэффициент адаптации на некоторой достаточно малой величине при , т.е. после выполнения шагов коэффициент адаптации должен оставаться постоянным. Величина выбирается таким образом, чтобы процесс «активной адаптации» можно было считать законченным.
Б.Уидроу [10], рассматривая процедуру адаптации, фиксирует момент, когда вектор параметров (в наших обозначениях - вектор //) фактически перестает меняться. Время активной перестройки параметров («время обучения») различно для разных временных рядов и, соответственно, для различных моделей АРСС. Однако можно, рассмотрев на моделях разных процессов процедуру адаптации, выбрать величину заведомо большую, нежели максимальный период обучения.
Таким образом, можно определить величину коэффициента адаптации как функцию такта s адаптации следующим образом:
1. начальное значение, т.е. р 0 положить равным
1;
2. на последующих тактах «обучение» модели
полагать р'=-, к = 2 ,. . . d;
к к '
3. после достижения такта d полагать р постоянным, равным -.
Процесс адаптации прогнозирующих моделей иллюстрируется следующим простым примером.
Был сгенерирован временной ряд по заданной модели авторегрессии первого порядка (АР1). При генерации был задан коэффициент 0! равный 0,8. Начальное значение ряда уг полагалось равным нулю. Прогноз на следующий такт вычислялся по формуле yfe = 0 - -yfe_ то есть для момента 2 прогноз равен нулю. К прогнозу добавлялся «шум» at, взятый из выборки, соответствующей нормальному распределению с нулевым математическим ожиданием и средним квадратическим отклонением <т (СКО) равным 0,1. Был сгенерирован ряд длиной 400 точек. Такая процедура соответствует известной схеме пропускания белого гауссового шума через формирующий фильтр. Фрагмент такого ряда представлен на рис. 2.
0.5
Y1i 0
-0.5
0 100 200 300 400
Рис. 2. Фрагмент временного ряда (процесс авторегрессии 1 рода)
Сгенерированный ряд экспортировался в ППП STATISTICA, где по данному ряду строилась его математическая модель. Естественно, что полученная модель является моделью авторегрессии первого порядка с коэффициентом ф1= 0,825. Будем считать это значение «истинной» величиной параметра авторегрессии. Средний квадрат остатков (т.е. оценка остаточной дисперсии) оказалась практически равной 0,01 (с точностью до 1 ■ 10-4), что соответствует СКО возбуждающего шума. Будем считать эти характеристики «ориентирами» для последующего процесса адаптации.
Затем аналогичная процедура была проделана для выборки из того же временного ряда (первые 50 точек). Получена оценка параметра АР (ф1=0,850). Далее включалась процедура адаптации, соответствующая формуле (13). Операция проектирования в этом случае сводится (при необходимости) к ограни-
чению по модулю |1| величины коэффициента авторегрессии.
На рис. 3 представлен процесс адаптации коэффициента 0 Длина временного ряда составляет 350 точек, т.к. первые 50 точек были использованы для получения начальной оценки ф1.
0.95 -
0.9
Xi
0.85
0.8 -
0 100 200 300 400
i
Рис. 3. Пример адаптации коэффициента ф,
Рис. 4. Начальный этап процесса адаптации
Из результатов видно, что установившееся значение оценки ф1 равно примерно 0,825. То есть можно
считать, что «обучение» модели прошло весьма успешно. Укрупненный фрагмент этого процесса (начальный этап адаптации) приведен на рис. 4. Из рисунков видно, что «обучение» можно считать законченным после 50 шагов.
Квадраты «остатков» от прогнозов представлены на рис. 5. Их средняя величина равна 0,01, что опять же соответствует среднему квадратическому отклонению возбуждающего шума. Следовательно, шум градиента при подобном выборе коэффициента р в стационарном режиме практически исключен.
i
Рис. 5. Квадраты «остатков» от прогнозов
Нестационарности, которые могут иметь место при работе с прогнозирующими моделями, следует рассматривать как разладку моделей. При этом крайне важно выбрать соответствующие критерии и при обнаружении разладки устанавливать коэффициенты р в «1».
Необходимо также решить задачу нахождения области допустимых значений параметров АР и СС при порядках АР и СС, больших двух (в литературе ограничиваются рассмотрением области допустимых значений только для моделей второго порядка).
Все эти вопросы требуют специального исследования и в данной работе не рассматриваются.
1. Браммер К., Зиффлинг Г. Фильтр Калмана-Бьюси / пер. В.Б. Колмановского. М.: Наука, 1982. 198 с.
2. Гантмахерт Ф. Р. Теория матриц. 2-е изд.доп. М.: Наука, 1966. 576 с.
3. Филимонов Б.П., Хрусталев Ю.П. Определение значений частот кванто-механических стандартов частоты адаптивными методами // Измерительная техника. 1987. №7, С. 27-30.
4. Эльясберг П. Е. Определение движения по результатам измерений. М.: Наука, 1976. 416 с.
5. Хрусталев Ю.П., Акулов В.М., Ипполитов А.А., Курыше-ва Л.Н., Обработка данных, полученных по результатам взаимных измерений вторичного эталона времени и частоты // Вестник Иркутского государственного технического университета. 2012. № 7. С. 22-29.
6. Ипполитов А.А, Хрусталев Ю.П. Субоптимальная филь-
ский список
трация в системах с неполной матрицей наблюдений // Информационные и математические технологии в науке и управлении: тр. XV Байкальской Всерос. конф. Иркутск: ИСЭМ СО РАН, 2010. Ч. I. С. 174 -182.
7. Бокс Дж. Анализ временных рядов. Прогноз и управление / пер. Д. Бокс, Г. Дженкинс. М., 1974. Вып.1. 406 с.
8. Ермольев Ю. М. Методы стохастического программирования. М.: Наука, 1976. 239 с.
9. Хрусталев Ю.П. Статическая и динамическая обработка данных, получаемых в процессе ведения эталонов времени и частоты // Измерительная техника. 2004. №6. С. 20-23.
10. Стационарные и нестационарные характеристики обучения адаптивных фильтров, использующих критерии минимума СКО / Б. Уидроу, Дж. М. Маккуи, М.Г. Ларимор, С.Р. Джонсон // ТИИЭР. 1976. Т. 64, №8. С. 37-51.
