Построение динамических стохастических моделей систем с неполной матрицей наблюдений Текст научной статьи по специальности «Математика»

At 2

Ц f (x,t)e 2 ■ cosQkt ■ Xk (x)xdt =

1V1

Г(( + gkvk )=1.

(16)

0 0

D

(13)

k=1

2k Q ■

k '

Задача (15)-(16) - задача на условный экстремум с соответствующей функцией Лагранжа. Поскольку

1 1

Ц f (x,t)e 2 ■ sin Qkt ■ Xk (x~)dxdt =

i

J sin2 knxdx

00

D

(14)

1

2'

1k

■Q

то в (15) минимизировать при условии (16) нужно одномерный интеграл

12

2 J

k=1

Обозначим правые части моментных равенств (13) и (14) соответственно ]лк и Ук.

В пространстве Ь2 (о), когда

||/(х,1 ) < I, I > 0, конечномерная проблема момен- ( ) .„ .

! - *(х^(^ ) тов для функции ] (х,1) , минимальной, например, по

М i

Г (dk ■cos Qkt+gk ■sin Qkt V

Akt

dt

Решение задачи (1)-(5) получаем в виде

M

(

k=1

dk cos Qkt +

V+ gk sin Qkt

Л

e 2 sin km

kl J

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

(

min

1 1 JJ

'dk ■ cos Qkt

Г

,gk ■sin Qkt

k j

k=1 Ail xe 2

■ Xk (x)

x sign

М AtL

Г/dk* cosQkt + gk* sinQkt)e 2 sinknx

k=1

Здесь dk*' gk*' k = 1,2'... - те числа, при которых

dxdt = Я

(15)

достигается условный минимум, равный

1

( Я )2

, вы-

при условии, что числа dk и gk удовлетворяют соот- ражения (15).

ношению

При р+ 2 и р > 1 задача (15)-(16) приводится к решению нелинейной системы уравнений.

Библиографический список

1. Бутковский А. Г. Методы управления системами с распре- 2. Мальцев Л.Е. Приближённое решение некоторых динами-делёнными параметрами. М.: Наука, 1975. 586 с. ческих задач вязкоупругости // Механика полимеров. 1978.

№ 2. С. 210-218.

T 1

2

At

k

2

X

УДК 681.3:529.7

ПОСТРОЕНИЕ ДИНАМИЧЕСКИХ СТОХАСТИЧЕСКИХ МОДЕЛЕЙ СИСТЕМ С НЕПОЛНОЙ МАТРИЦЕЙ НАБЛЮДЕНИЙ

Ю.П. Хрусталев1

Национальный исследовательский Иркутский государственный технический университет, 664074, г. Иркутск, ул. Лермонтова, 83.

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

Ключевые слова: динамические стохастические модели; метод наименьших квадратов; функционал; недооп-ределённая система; групповые эталоны.

1Хрусталев Юрий Петрович, кандидат технических наук, доцент кафедры вычислительной техники, тел.: (3952)405107, (3952)405163, (3952)510698.

Khrustalev Yury, Candidate of technical sciences, Associate professor of the chair of Computer Science, tel.: (3952) 405107, (3952) 405163, (3952) 510698.

CONSTRUCTION OF DYNAMIC STOCHASTIC MODELS FOR SYSTEMS WITH INCOMPLETE OBSERVATION MATRIX

Yu.P. Khrustalev

National Research Irkutsk State Technical University, 83, Lermontov St., Irkutsk, 664074.

The author deals with the questions of construction of dynamic stochastic models for systems with incomplete observation matrix. It is shown that parameter estimation of autoregression models - integrated moving average- can be obtained by minimizing the sum of the squares of measurement results divergence from their projections, calculated on the basis of these models. The author proposes the procedure to construct dynamic stochastic models in the presence of linear trends in time series. 8 sources.

Key words: dynamic stochastic models; method of least squares; functional; underdetermined system; group standards.

Динамическими называются системы, состояние которых зависит от времени. Будем считать, что /-е состояние описывается параметром у. Тогда состояние многомерных динамических систем может быть представлено ^-мерным вектором состояния 4(7). В задачах оценивания состояния линейных динамических систем наблюдаемая векторная величина г^) связана с вектором матрицей наблюдений А^):

гю = А(1) ■ У() (1)

Размерности векторов г и У в общем случае не совпадают. В реальных технических системах матрица А, как правило, не зависит от времени, т. е. А(^=А. Если размерность вектора г больше размерности А, то имеет место недоопределённая система, в которой число измерений больше числа оцениваемых параметров. Для получения оценок вектора состояния У в момент t в этом случае используются статистические методы, сводящиеся к минимизации суммы квадратов

«невязок» между вектором гА) и его оценками 2(7) [5].

В недоопределённых системах размерность выбора г меньше размерности вектора У, данная ситуация имеет место в групповых эталонах, например, в эталонах времени и частоты [8]. Для получения оценок вектора состояния У в момент ^ необходимо привлечение некоторой априорной информации о состоянии системы, либо введение дополнительных ограничений на вектор У (например, в качестве решения

уравнения (1) берется вектор 4 с минимальной нормой [2]). Тогда оценка 4 может быть найдена с помощью псевдообратной матрицы А+ :

/ = А+ • г . (2)

Для динамических систем априорная информация чаще всего представляется в виде прогнозов, вычисленных для моментов ^ на основе законов «движения» системы, т. е. её динамики. В технических системах, таких как эталоны времени и частоты, вектор состояния обычно оценивается через равные промежутки времени.

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

измерения, выполняемых в эталонах времени и частоты, описана в [6]. В настоящее время существуют пакеты прикладных программ для построения моделей АРПСС, в которых реализована методика Бокса -Дженкинса [1].

Особенность рассматриваемой задачи состоит в том, что в распоряжении исследователя отсутствуют временные ряды 4(7), для которых нужно построить математические модели (имеются только результаты наблюдений 2,(1)). Поэтому предложено для построения моделей АРПСС, описывающих не составляющие вектора 4(7) -уь где i = 1..п , минимизировать функционал:

j=1

(3)

где N - длина временного ряда; 2) - прогноз вектора измерений для момента tj.

В групповых эталонах результаты измерений представляют собой разности между значением величины, воспроизводимой некоторым из элементов эталона (опорный элемент) и другими элементами, т. е. 2(7) = уоп — yi. Пусть в качестве опорного выбран первый элемент, тогда вектор наблюдений определён

как

2T(t) = [Уг - У2.У1 - Уз>-'У1 - Уп1 Матрица измерений в этом случае имеет вид:

(4)

A =

1 -10 10 -1 100

0 0 -1

Нетрудно увидеть, что при включении в групповой эталон п элементов, ранг матрицы А равен п — 1.

При субоптимальной фильтрации оценка состояния опорного элемента может быть найдена по формуле [8]:

Ул*)=i & (У г (1)+^ (t»,

(5)

где g, - вес /-го измерения;

У (1) -

прогноз состояния

/'-го элемента, вычисляемый на предыдущем такте обработки данных (т.е. прогноз на 1 шаг вперед).

1=1

При выполнении условия ^gi = 1 оценка

1=1

У1(1) - несмещённая.

Прогнозы Ус (1) вычисляются на основе моделей АРПСС:

У (1)=вв ■ х, (6)

где вв =[ф1, фр; вь-, ва] - вектор параметров

модели АРПСС ¡-го элемента; Xе =[ум, Ум,..., у-р;

-ам,...,-а-а] - вектор состояния прогнозируемого процесса для ¡-го элемента; р, ц - порядки авторегрессии и скользящего среднего для модели ¡-го элемента; ф -обобщенный коэффициент авторегрессии; а - ошибка ¡-го прогноза.

С учётом введённых обозначений функционал (3) принимает вид:

N п—1 ) )

—смо—Л>+1(0)]2- (7)

]=11=1

Поскольку прогнозы У (/) вычисляются на основе моделей АРПСС, функционал 3 зависит от векторов параметров моделей въ . Минимизация функционала

приводит к получе-нию оптимальных оценок параметров моделей АРПСС для каждого из элементов группового эталона - в*.

Таким образом, алгоритм нахождения оптимальных значений параметров моделей АРПСС заключается в следующем:

1. Ввод размерности решаемой задачи, т.е. величин п, Ы, б, р, ц (б - порядок разности; р - порядок АР; ц - порядок СС).

2. Ввод результатов измерений , начальных значений векторов вО.

3. Ввод начальных оценок У

4. Вычисление оптимальных значений векторов

ß'r

a) вычисление функционала 3 (ß).

b) вычисление градиента V3(ßk ) в точке

(ßk) (k=0,1,2...). Если норма градиента меньше заданной величины, определяющей точность решения оптимизационной задачи - решение найдено. При

этом оценка i-ой составляющей вектора состояния (Y i) находится из выражения

Y=Y -, (8)

где Y - оценка состояния опорного элемента, найденного по формуле (5);

c) решение задачи одномерной оптимизации (поиск значения /, удовлетворяющего условию

3(Ä,ß,g) = min 3[(ß,g) + ÄVk(ß,g)J);

с1) нормировка весов д (в начале оптимизационной процедуры веса д, (формула (5)) полагаются равными,

о 1

так что qi = —. Затем, после нахождения текущих п

оценок векторов @к, на к-ом шаге веса выбираются

обратно пропорциональными остаточным дисперсиям или средним квадратом ошибки прогнозов для каждого временного ряда У, );

е) переход на этап 4,а.

Работоспособность алгоритма при известной структуре моделей АРПСС (порядки р, б и ц) и отсутствии линейных трендов (и трендов более высокого порядка) во временных рядах продемонстрирована в работе [7].

Методы идентификации структуры моделей АРПСС по имеющимся временным рядам были предложены Боксом и Дженкинсом [1]. Эти учёные на основании анализа большого количества различных временных рядов установили, что порядки авторегрессии (р) и скользящего среднего (ц) для рядов, отражающих функционирование реальных систем, не превышает 5, а в большинстве случаев равны 2 или 3.

При разработке алгоритмов оценивания относительных отклонений частоты водородных стандартов [5] были построены модели указанных временных рядов. При этом выводы Бокса и Дженкинса подтверждались. Следует учитывать то обстоятельство, что, в отличие от задачи, решаемой в рамках данной работы, ряды относительных отклонений частоты водородных стандартов считались известными. Эти данные были получены на основе метрологического анализа результатов наблюдений, в процессе которого привлекалась дополнительная информация о состоянии эталона, получаемая в процессе «внешних сличений» с другими эталонами времени и частоты, а также использовались различные неформализованные процедуры, основанные на опыте высококвалифицированных специалистов-метрологов.

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

1. Рассматриваются «простейшие» (тривиальные) модели, в которых порядки р и д равны нулю.

2. Последовательно, начиная с опорного элемента, увеличиваются на единицу порядки авторегрессии, так что они возрастают с 0 до 3-х. Для каждого случая «прогоняется» описанный выше алгоритм и фиксирует

ются значения функционала .

3. Аналогичная процедура проводится и для порядков скользящего среднего д.

4. Выбираются в качестве конечных результаты, при которых «средний квадрат» (т.е. значения функционала, делённые на число оцениваемых параметров) минимален.

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

поскольку независимость остатка (а,- = z(г,)-))

у ' у

неочевидна, так как в каждом результате наблюдений содержится информация о состоянии опорного элемента. Следует также иметь в виду, что при незначительной длине временных рядов (N«100) специализированный программный комплекс работает весьма быстро, и затраты машинного времени крайне незначительны. Поэтому можно использовать рассмотренную выше процедуру простого перебора величин р и Я-

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

• введением в состав вектора параметров моделей АРПСС свободных членов ва ;

• включением в модель векторной детерминированной функции ф(1).

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

Рассмотрим наиболее простые тренды - линейные, которые весьма часто встречаются в технических системах.

Если при минимизации функционала (7) рассматривать расширенный вектор параметров моделей АРПСС

Рв =[ср1,ф2,...,фр,е1,е2,...,еч\ео ],

где ео - свободный член, который при порядке разности d = 1 описывает детерминированный линейный тренд, то оценки векторов Рр не всегда сходятся к своим истинным значениям. Этот результат становится очевидным, если учесть, что при неполной матрице наблюдений мы имеем дело с недоопреде-лённой системой. Действительно, в этом случае можно найти бесчисленное множество решений линейных алгебраических уравнений. Поскольку в моделях АРПСС коэффициенты взаимосвязаны (увеличение одних коэффициентов можно компенсировать уменьшением других [5]), то это обстоятельство и объясняет случаи отсутствия сходимости процесса построения модели.

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

1. На первом этапе попытаемся найти оценки параметров линейных трендов.

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

3. Оценки параметров моделей авторегрессии -скользящего среднего находятся минимизацией функционала °.

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

При рассмотренной выше матрице наблюдений А (размерность матрицы (п — 1) х п) и линейных трендах вида

У г = Ь0, + Ь1, ■ Г . (9)

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

zi = (Ъо1 — Ъ0с) + (Ьи — Ьн) ■ Г, (10)

где г = 1,2,3,...,п; Г = 0,1,2,..., N.

Коэффициенты Ь01 находятся из начальных условий задачи (предполагается, что в момент 1=0 эти величины известны из результатов «внешних сличений» группового эталона). Таким образом, проблема состоит только в нахождении коэффициентов Ь,, г = 1,2,3,...,п . Если будет найдено хотя бы одно какое-либо значение коэффициента Ь1,, остальные,

при известных коэффициентах Ь01, немедленно находятся из системы уравнений (10).

Если результаты наблюдений

г = 0,1,2,...^.:

г = 1,2,3,...,п не содержат детерминированной линейной составляющей, то это может означать, что линейные тренды отсутствуют во временных рядах, соответствующих опорному элементу и элементу с номером / (возможностью существования в этих временных рядах линейных трендов с одинаковыми углами наклона в реальных системах можно пренебречь ввиду малой вероятности такого события). Таким образом, задачу вычисления оценок параметров детерминированных линейных функций можно сформулировать как задачу проверки гипотезы Н0: коэффициент Ь1 линейной функции, содержащейся в каком-либо временном ряде г,, равен нулю.

Если данная гипотеза не отвергается хотя бы для одного временного ряда то можно считать, что

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

При наличии линейного тренда во временном ряде, представляющем состояние опорного элемента, гипотеза Н0 будет отвергнута для всех временных рядов 2^ . В этом случае необходимо поменять местами

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

р = 1 с„2 ■в'(1 — в)"

(11)

1 = 2

где п - число элементов в групповом эталоне; сп -

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

Очевидно, что с ростом п растёт вероятность обнаружения линейных функций yi(t). Например, при одинаковой вероятности наличия и отсутствия линейных трендов (т.е. при в =0,5) вероятности их идентификации (как функции п) равны: п =2; р=0,5; п =5; р=0,78; п=7; р=0,94.

Таким образом, с ростом числа элементов, входящих в групповой эталон, растёт не только его техническая надёжность, но и «алгоритмическая надёжность».

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

следует найти временной ряд 2^1), для которого

вероятность наличия тренда минимальна. Это сделать не трудно, так как подобная задача сводится к нахождению вероятностей распределения Стьюдента с п-1 степенью свободы по вычисленным значениям величины 1

N 1

)2 72

(12)

^ = ъ1[ I (1 — )212/$, 1=1

где Ъ = I(1 — ) ■ (2, — 2)/ I (/ — 1р )2 ;(13)

1=1

1=1

1ср = (11)/2;

2 = 11 2 г

N

(14)

(15)

N ) 2

5 = (I (2. — 2 ))/(N — 2) ;

(16)

1=1

2 = Ъ0 + Ъ ; 1 = 1,2.

(17)

Соотношения (12)-(17) соответствуют формулам парной регрессии [7], когда независимая переменная принимает значения 1 = 0,2,...,N, а откликом являются величины г¡.

Если из каких-либо соображений известны начальные значения углов наклона линейных трендов, присутствующих во временных рядах (например, по результатам обработки данных на предыдущем интервале, включающем N тактов), то можно вычислить начальные оценки параметров детерминированных функций, входящих во временные ряды Последующая процедура заключается в проверке гипотезы: коэффициент Ь, линейной функции, содержащейся в каком-либо временном ряде, равен его начальному значению.

Если хотя бы для одного временного ряда эта гипотеза не отвергается, можно полагать, что углы наклона линейного тренда у опорного элемента и элемента с номером , остались без изменений. Учёт этого обстоятельства делает систему уравнений (1) разрешимой (если речь идет об идентификации линейных функций). «Двухшаговая» процедура определения коэффициентов Ь¡, определяющих углы наклона линейных функций, содержат погрешности двух видов:

• погрешности оценивания углов наклона рядов измерений и связанные с ними вероятности ошибки при проверке гипотез об отсутствии линейных трендов в исходных временных рядах;

• погрешности, вызванные заменой минимального угла наклона (в случае, когда гипотеза об отсутствии линейных трендов отвергается) нулевым значением [формулы (12)-(17)].

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

Рассмотрим функционал 3 :

N п—1 ) ) 2

3 = 11[ 2Р—(УЦ — У 1.1+1)].

1=1 1=1

Пусть, для наглядности, имеется всего два исходных ряда у. Следовательно, в нашем распоряжении только один ряд измерений. Пусть также ряды у, (¡=1,2) не содержат случайной составляющей. Тогда

2р = Ур — У]2. Прогнозы для момента у будут числяться по формулам:

вы-

1=1

Ai = ¿oi + bii ■j;

^j 2 = b02 + b12 ■j; j = 2,3,.. N.

В групповых эталонах начальные значения состояний b0¡ элементов известны из результатов внешних сличений. Поэтому коэффициенты b0í можно полагать равными нулю, с последующими введениями поправок в результаты вычислений. Для простоты

изложения обозначим в11 как /3¡, а в12 как в2, а

их оценки, соответственно, как в11 и в2. Тогда функционал (7) примет вид

3 =

j=

Оптимальные

дЗ „ дЗ

S[ (ß ß2 ) •j -(ß -ß) •

j

ß = 0

ß

и

dß2

оценки найдутся из условии = 0, что приводит к системе линей-

ных уравнении:

N

Si ßi-ßi) - (h - b2)\( -j2) = 0;

j=i

N

Si ßi-ß?) - (bi - b2)]-j2 = 0,

(18)

j=i

(в (18), в соответствии с общепринятым в МНК подходе, Д. заменены на (/=1,2)).

Оценки Ьг можно найти, вынося выражения, не

связанные с ], за знак суммы. Тогда условия оптимальности примут вид:

-(ßi-ß2) + (bi -b2) = 0; (ßi-ß2) - (bi - b3) = 0.

(19)

Система уравнений (19) - линейно зависимая, следовательно, не имеет единственного решения, чего и нужно было ожидать. Поэтому при использовании функционала 3 для поиска оптимальных оценок линейных трендов крайне важно выбирать начальные

оценки , для чего и рекомендуется использовать

описанную выше «двухшаговую процедуру».

Проверка работоспособности процедуры минимизации функционала (7) при отсутствии детерминированных функций во временных рядах проводилась методом моделирования и дала положительные результаты. Для поиска оптимальных оценок Д. можно

использовать тот же самый программный комплекс, заменив прогнозирующую функцию (вместо моделей АРСС вводить линейные модели).

Библиографический список

1. Бокс Д., Дженкинс Г. Анализ временных рядов. Прогноз и управление /пер. с англ. под ред. В.Ф. Писаренко. М.: Мир, 1974. Кн. 1. 406 с.

2. Гантмахер Ф.Р. Теория матриц: учеб. пособие для вузов. М.: Наука,1966. 576 с.

3. Дрейпер Н., Смит Г. Прикладной регрессионный анализ: в 2-х кн., 2-е изд., перераб. и доп. / пер. с англ. Кн. 1. М.: Финансы и статистика, 1986. 365 с.

4. Кашьяп Р.П., Рао Р.П. Построение динамических стохастических моделей по экспериментальным данным. М.: Наука, 1985. 384 с.

5. Лоусон И., Хенсон Р. Численное решение задач методом наименьших квадратов. М.: Наука, 1986. 232 с.

6. Хрусталев Ю.П., Спиридова Е.В. Алгоритмы обработки измерительной информации, получаемой в процессе хранения единиц времени и частоты // Техника средств связи. Серия «Радиотехнические измерения». М., 1986. Вып. 1. С. 58-72.

7. Хрусталев Ю.П., Овечкина А.А., Щербаков Е.В. Построение моделей временных рядов по результатам наблюдений в динамических системах. В сб. Методы исследования и моделирования технических, социальных и природных системах. Новосибирск: Наука, 2003. С. 293-307.

8. Хрусталев Ю.П. Статистическая и динамическая обработка данных, получаемых в процессе ведения эталонов времени и частоты // Измерительная техника. 2004. № 6. С. 20-24.