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

2011

ВЕСТНИК ТОМСКОГО ГОСУДАРСТВЕННОГО УНИВЕРСИТЕТА Управление, вычислительная техника и информатика

№ 4(17)

ОБРАБОТКА ИНФОРМАЦИИ

УДК 519.21

А.М. Горцев, В. Л. Зуевич

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

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

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

Важной сферой приложения теории массового обслуживания является проектирование и создание информационно-вычислительных сетей и различных сетей связи, которые можно объединить единым термином - цифровые сети интегрального обслуживания (ЦСИО). Данная сфера была определена развитием информационных технологий в конце ХХ века. Возникла необходимость в разработке математических моделей потоков событий, адекватно описывающих реальные информационные потоки, функционирующие в ЦСИО. Одними из первых работ в этом направлении были [1-3]. Отметим, что на практике параметры, характеризующие поток событий, частично либо полностью неизвестны. Кроме того, параметры могут изменяться с течением времени случайным образом, что приводит к рассмотрению дважды стохастических потоков событий. Поскольку функционирование системы обслуживания непосредственно зависит от параметров входящего потока, важной задачей является оценка в произвольный момент времени его параметров по наблюдениям за этим потоком. Исследования по оценке параметров дважды стохастических потоков были проведены, например, в работах [4 - 6].

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

1. Постановка задачи

Рассматривается асинхронный дважды стохастический поток событий с произвольным конечным числом состояний (далее асинхронный поток либо просто поток). Интенсивность потока является кусочно-постоянным случайным процес-

сом Х(/) с п состояниями: Хь Х2, ... ,Хп (Х1 > Х2 > ... > Х„ > 0). Процесс (поток) в момент времени / находится в г-м состоянии, если Х(/) = Хг (г = 1, п). В течение времени пребывания в г-м состоянии поток ведет себя как пуассоновский с интенсивностью Хг (г = 1,п ). Длительность пребывания в г-м состоянии есть экспоненциально распределенная случайная величина с функцией распределения

П ___ ___

(т) = 1 -га‘‘т, где агг = - ^ аг]- (г = 1,п ); аг] > 0 (г,] = 1,п , г ^ ]) - интенсив-

]]*г

ность перехода процесса Х(/) из состояния г в состояние ], т.е. величины а,]- образуют матрицу интенсивностей (матрицу инфинитезимальных коэффициентов) переходов между состояниями ||аг] Щ . В сделанных предпосылках Х(/) - транзитивный марковский процесс [7].

Значения параметров потока Хг, аг] (г, ] = 1, п , г ^ ]) неизвестны. Процесс Х(/) является принципиально ненаблюдаемым. Предполагается, что о потоке известно только число состояний п и наблюдению доступны только моменты наступления событий потока /ь /2, ... . Необходимо по наблюдениям /ь /2, ... оценить параметры потока Хг, а] (г, ] = 1, п , г ^]) в момент окончания наблюдения за потоком.

Рассматривается стационарный (установившийся) режим функционирования наблюдаемого потока событий, поэтому переходными процессами на интервале наблюдения (/0 , /), где /0 - начало наблюдений, / - окончание наблюдений (момент вынесения решения), пренебрегаем. Без потери общности можно положить /0 = 0. Обозначим 0 = (Х1,...,Хп,аг]-;г,] = 1,п,г * ]) - вектор неизвестных параметров

потока, 0(/) - вектор соответствующих оценок параметров в момент времени /. 0к и 0к (/) - к-е компоненты вектора параметров и вектора оценок соответственно. Обозначим N - размерность векторов 0 и 0(/), N = п2. Обозначим р( 0 | /) = = р( 0 | /ь /2, ..., /т) - апостериорная плотность распределения вероятностей вектора параметров 0 в момент времени / при условии, что в моменты времени ^, ^ ..., 4; (0 < /1 < /2 < ... < /т < /) наблюдались события потока. Обозначим © = {Х1 > Х2 > . > Хп > 0, аг] > 0; г,] = 1,п , г' ^]} - область значений вектора параметров 0 . Будем

использовать оценку 0(/), оптимальную в смысле минимума среднеквадратического отклонения от истинного значения вектора 0 :

0к (/) = 10кр(01 /)й0 (к = ),

©

где й0 = d0\d02 ... d0N, или в векторной форме

0(/) = 10р(0 | /)й0. (1)

©

Выражение (1) дает оптимальную оценку параметров потока в виде апостериорного среднего. Чтобы воспользоваться формулой (1), необходимо найти выражение для апостериорной плотности р( 0 | /).

2. Явный вид апостериорных вероятностей состояний потока

Для нахождения апостериорной плотности р( 0 | /), как будет показано в разделе 3 настоящей статьи, необходимо знать ю(Ху | /) (у = 1, п) - апостериорные вероятности пребывания асинхронного потока ву-м состоянии в произвольный момент времени t. Явный вид апостериорных вероятностей ю(Ху | /) (у = 0, п) для любого момента времени найден в работе [7].

Согласно [7], поведение ю(Ху | 0 (у = 1, п) на полуинтервале [4 , tk+1) (к = 1, 2, ...) между соседними наблюдёнными событиями асинхронного потока и на полуинтервале [4, и) между началом наблюдений и наблюдением первого события определяется выражением

«кУ^(Мк) _

“У) =-ПТ-п---------------- (У = 1, п, к = 0, 1, ... ), (2)

¿¿зд (/к )ею(‘-к)

- =1 I=1

пп

где 4 < ^ < 4+ь ю(Ху | tk) = ю(Х; | tk + 0) = ¿(1к^ г, (4) = ¿.й~1ю(Хг | tk + 0),

г=1 1=1

к = 0, 1, ...; ю, (, = 1,п) - собственные числа матрицы ^ = |ЦД”, й, = а*- (г ^ /),

йп = а,, - X, (г,, = 1,п ); я, - элементы матрицы £ = ||яг7|Щ , в которой ,-й столбец является собственным вектором, соответствующим собственному числу ю, ; я-1 (I, г = 1, п ) - элементы матрицы £-, обратной матрице £. В момент времени (к (в момент наступления события асинхронного потока) апостериорные вероятности (2) претерпевают разрывы 1-го рода [7] и пересчитываются по формуле

Х ую(Х у | Гк - 0) —

ю(Ху\к + 0) =-^----}—к---- (у = 1, п, к = 1, 2, ... ), (3)

¿хгю(х,ик - 0)

г=1

где ю(Ху | ^ - 0) вычисляется по формуле (2) в момент t = 4, когда t изменяется в полуинтервале [4ч , 4), соседнем с полуинтервалом [4 , 4+1). В качестве начального значения ю(Ху | t0 + 0) = ю(Ху | t0 = 0) в (2) выбираются априорные финальные вероятности Пу (у = 1, п) состояний процесса Х(^ в стационарном режиме ( ^ да), которые удовлетворяют системе линейных алгебраических уравнений [7]

¿Пгау = 0 ( 1 = 1п - 1 ), ¿Пу = 1. ()

г=1 1=1

3. Нахождение явного вида апостериорной плотности вероятностей и оценок параметров потока

Получим рекуррентную формулу для определения апостериорной плотности вероятностей. Воспользуемся методикой [8]: рассмотрим наблюдения за потоком через равные достаточно малые промежутки времени длительности At, а затем устремим At к нулю. Пусть наблюдения за потоком начинаются в момент времени

ґ = 0 и время ґ изменяется дискретно с конечным шагом Дґ: ґ = кДґ (к = 0,1,...). Обозначим г(кДґ) - число событий, наблюдённых на полуинтервале времени [(к - 1)Дґ, кДґ) (г(кДґ) = 0,1,., к = 0,1,.). На полуинтервале [-Дґ, 0) наблюдений не производится, поэтому г(0) можно положить произвольным, например, положим г(0) = 0. Обозначим г (тДґ) = (г (0), г (Дґ),..., г (тДґ)) - последовательность значений количества наблюденных событий на временных полуинтервалах [(к - 1)Дґ, кДґ) (к = 0, т). Рассмотрим момент времени ґ, такой, что ґ = тДґ, ґ + Дґ = = (т + 1)Дґ. Тогда имеем г(тДґ) = г(ґ), г((т + 1)Дґ) = г(ґ + Дґ), г (тДґ) = г (ґ), Г ((т + 1)Дґ) = г(ґ + Дґ).

Введем р( 0 | г (тДґ)) = р( 0 | г (ґ)) = р( 0 | ґ) - апостериорную плотность вероятностей вектора параметров 0 в момент времени ґ при условии, что на полуинтервалах времени [(к - 1)Дґ, кДґ) (к = 0, т) наблюдалось г(кДґ) событий потока соответственно. Аналогично введем р( 0 | ґ + Дґ) - апостериорную плотность в момент времени ґ + Дґ.

Теорема 1. Апостериорная плотность р( 0 | ґ + Дґ) определяется рекуррентной формулой

где ю(Х,-1 ґ) (у = 1,п) - апостериорная вероятность того, что поток в момент времени ґ (ґ > ґ0) находится ву-м состоянии, определяемая формулами (2), (3). Доказательство. Используя формулу для условной вероятности, находим

Используя в последней дроби формулу полной вероятности для условной вероятности, получим

Рассмотрим в (6) произведение р(г^ + А^| г^),0,X(t) = Ху)р(Х(0 = Ху1 г(0,0). Имеем p(X(t) = Ху | г(t),0) = ю(Ху 11) (у = 1,п), где ю(Ху- | t) (у = 1,п) - апостериорная вероятность того, что поток в момент времени t находится в у-м состоянии, определяемая в зависимости от момента времени t формулами (2), (3). В силу предпосылок, количество наблюдённых на полуинтервале [, t + А^ событий не зависит от последовательности г ^) наблюдённых до момента t событий, а зависит только от состояния потока в момент времени t, т.е. от значения интенсивности наступления событий потока Х(^ = Ху (у = 1, п) в момент времени t. Поэтому имеем р(г^ + А^| г^),0,Х^) = Ху) = р(г^ + А^| 0,X(t) = Ху). Поскольку асин-

г (ґ+Дґ )

р(01 ґ + Дґ) = р(0 | ґ)

(5)

р(01 г (ґ + Дґ))

р(0, г(ґ + Дґ), г (ґ)) _ р(г(ґ + Дґ) | 0, г (ґ))р(0 | г (ґ))р(г (ґ))

р(г (ґ + Дґ))

р(г (ґ + Дґ))

р(01 г (ґ+Дґ)) = р(01 г (ґ ))—-

р(г (ґ)) ^

р(г(ґ+Дґ))

р(г(ґ+Дґ)|г(ґ),0Д(ґ)=Х;)р(Х(ґ)=Ху |г(ґ),0) .(6)

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

_ (X At)Г(Î+AÎ) .д .

p(r(t + At)| 9,X(t) = X.) = —--e 1 . Таким образом, находим формулу (6)

1 r (t + At )!

в виде

p(9 1 r(t + At)) = p(9 1 r(t)) p(r +A ю(Х 1 11)( J ^ e~XjAt.

p(r (t + At)) 1=1 r (t + At)!

TT p(r:(t ))

Находя в последнем выражении множитель ------------------------ из условия нормировки

Р p(r (t + At )) J P P

| p(9 | r(t + At))d9 = 1 , приходим к (5). Теорема доказана.

0

Рассмотрим случай, когда на полуинтервале [t, t + At) нет событий, т.е. r(t + At) = 0. Это означает, что полуинтервал [t, t + At) находится на временной оси между моментами наступления соседних событий tk и tk+1 (к = 1, 2, ...) либо между моментом начала наблюдений за потоком t0 и моментом наступления пер-

вого события Введем обозначение s(t, 9) = Xх 1ю(х 111).

1=1

Лемма 1. Апостериорная плотность р( 9 | /) между моментами наступления событий удовлетворяет интегро-дифференциальному уравнению

ф(91 /)

dt

■ = - p (9 11 )

s(t, 9) -1 s (t, 9) p (911 )d 9

_ 0

Доказательство. Поскольку r(t + At) = 0, получаем (5) в виде

(tk < t< tk+i, к = 0, 1, ...). (7)

X“(Xj 11)e

р(91 t+д/) = р(91 ¿)----]-=—п----------------

| р(9 | о£ю(Х 11 Д^9

© 1=1

или, разложив экспоненты в числителе и знаменателе в ряд с точностью до о(Д/), получим

1 - At| s(t, 9)p(9 11)d9 + o(At)

р(911 + Д[) = _ р(911) - s(t, 9) р(911 )Дt + о^)]

_ ©

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

р(9 11 + Дt) = р(9 11) -Дts(t, 9)р(9 11) + Дtp(9 11)| s(t, 9)р(9 11)ё9 + o(Дt).

0

Перенося в последнем равенстве р( 9 | ^ влево, деля на Дt и устремляя Дt к нулю, приходим к (7). Лемма доказана.

Асинхронный поток обладает свойством ординарности, поскольку в каждом состоянии ведет себя как пуассоновский. Поэтому вероятность наступления на полуинтервале [, t + Дt) двух и более событий равна о(Д^. Рассмотрим случай, когда на полуинтервале [X t + Дt) наступает одно событие потока в момент времени ^ ^ ^ < t + ДО. При этом г^ + Дt) = 1.

Лемма 2. Апостериорная плотность p( 9 | /) в момент наступления события пересчитывается по формуле

р(9 и, - 0)5(/, - о, 9)

р(0\/к + 0) =

| р(0 \ ґк - 0)і'(їк - 0,0)й9

(к = 1, 2, ...).

(8)

Доказательство. Поскольку г(/ + Д/) = 1, получаем (5) в следующем виде:

р(9 \ / + Д/) = р(9 \ /)-

і=1

Д/| р(9 | /)£ю(Я,] 11)Х^ я 1 9

0 1=1

или, разложив экспоненты в числителе и знаменателе в ряд с точностью до о(Д/) и вводя величины Д^ и Дt" такие, что t = ^- Д^, t + Дt = ^ + Д/", находим

Д/р(91 /к - д/')5(/к - д/', 9)+о(Д/)

р(0 \ Ч + Д/") =

Д/1 р(0 \ /к -Д/')5(/к -Д/', 0)й9+о(Д/)

Поделим числитель и знаменатель последней дроби на Д/, после чего устремим Д/ к нулю (при этом / = ґк - Д/ стремится к ґк слева, / + Д/ = ґк + Д/" стремится к ґк справа). В результате предельного перехода приходим к (8). Лемма доказана.

Замечание к лемме 2. В момент времени /0 начала наблюдений за потоком апостериорная плотность р( 0 \ /0) задается, исходя из априорных данных о параметрах асинхронного потока. Если таких данных нет, можно задать плотность р( 0 \ /0) как произведение N плотностей равномерно распределенных случайных величин, каждая из которых распределена в некотором интервале допустимых значений для соответствующего параметра потока.

Теорема 2. Поведение апостериорной плотности р( 0 \ /) на временной оси определяется интегро-дифференциальным уравнением (7) и формулой пересчета вероятностей (8), в которых /к < / < /к+1, р( 0 \ /к) = р( 0 \ /к + 0) (к = 0,1,.). В момент времени /0 апостериорная плотность р( 0 \ /0) задается согласно замечанию к лемме 2.

Доказательство следует из лемм 1 и 2 путем синхронизации формул (7) и (8).

Теорема доказана.

Теорема 3. Апостериорная плотность вероятностей р( 0 \ /) вектора параметров

0 на полуинтервале времени [/к, /к+1) (к = 0, 1, .) определяется формулой

р(0 \/к )ехР

р(0 \ /) = -

s(т, 0)й т

| р(0 \ ік )ехр

s(т, 0 )й т

(/к < / < /к+1, к = 0, 1, ...),

(9)

где р( 9 | /к) = р( 9 | /к + 0) вычисляется в момент наступления события /к (к = 1, 2, ...) по формуле (8), а в момент /0 апостериорная плотность р(9 | /0) задается согласно замечанию к лемме 2.

к

Доказательство. Преобразуем интегро-дифференциальное уравнение (7) к следующему виду:

ёр(91 ґ)

р(9 І ґ)

s(t, 9) -1 s(t, 9) р(91 ґ )ё 9

ж.

Интегрируя последнее интегро-дифференциальное уравнение в пределах от ґк до ґ, получаем

р(91 ґ) ґ __Г

_ р(9 | ґк) _ = і (к

1п

или, проделывая необходимые преобразования, получаем

і

| 5(т, 9) -| 5(т, 9)р(9 | т)ё9

Р(9 | ґ) = р(9 | ґк )ехр

5 (т, 9)ё т

ехр

115(т, 9)р(9 | т)ё9ёт

ґк ©

Находя последний сомножитель из условия нормировки | р(9 | ()й9 = 1, приходим

0

к (9). Теорема доказана.

Подставляя (9) в (1), получаем явный вид оценок параметров потока для произвольного момента времени наблюдения за потоком. При этом в момент времени 4+1 (к = 0, 1, ...) имеем

9«к+1) = 9(Ч+1 + 0) = |9р(9 | Гк+1 + 0)^9 , (10)

0

где р( 9 | 1к + 0) вычисляется по формуле (8).

4. Приближенные формулы для расчета оценок параметров потока

Сделаем важное замечание. При расчете вектора оценок параметров 9(/) по формулам (1), (8) - (10) возникают существенные сложности реализации алгоритма данного расчета на ЭВМ. Во-первых, интегрирование во всех указанных формулах ведется по Ж-мерной области 0, которая в общем случае может быть не ограничена. Во-вторых, интегралы являются Ж-кратными. Поэтому в настоящей

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

Идея приближенного алгоритма состоит в предположении достаточной близости вектора оценок 9(/) к истинному значению вектора параметров 9 , что позволяет воспользоваться разложениями в ряд Тейлора и значительно упростить формулы (1), (8) - (10) для вычислений.

В предположении достаточной близости оценок 9(/) к начальному значению

9(ґк + 0) (к = 0, 1, ...) разложим в интегралах вида|/(-,9)р(9 |ґ)ё9 .

входящих в

состав формул (1), (8) - (10), подынтегральные функции в ряд Тейлора в окрестности точки 9(ґк + 0) (к = 0, 1, ...), ограничиваясь первыми тремя членами разложения.

0

Введем обозначения:

т1 (г) = 9г (г) = 19гр(9 | г)й0 , с,, (г) = | (9, - т, ('))Ф1 - т1 (г))р(9 | г)й9,

/ (г, (к, 9) = ехр

,$(т, 9 )й т

/ ((, (к, 9) = 9,- ехр

,^(т, 9 )й т

(11)

где 9г, 9г (() - элементы векторов 9 , 9((), соответственно.

В (11) величины тг (() - апостериорные средние параметров 9г, г = 1, N. Величины с,1(') (,, I = 1, N) - ковариации компонент вектора параметров 9. Вычисление интегралов в формулах (11) представляет собой сложную задачу в части построения численного алгоритма на ЭВМ интегрирования по области 0.

Формула (1), согласно (9), с учетом обозначений (11), запишется в следующем виде:

р(91 % + 0)ехр

, (()=Р9<

5(т,9)^ т

I р(91 гк + 0)ехр

5(т,9 )й т

-а 9=

| /, (г ,гк, 9) р(91 гк + 0)й 9

0___________________________________

1/ (г,(к ,9) р (91 гк + 0)а 9

г = 1, N. (12)

Разложим в (12) функции (11) в ряд Тейлора в окрестности точки е('k + 0) = = т(гк + 0), пренебрегая членами порядка малости о(||л9||2), Д0 = 9-т(гк + 0). Получим

- Ж 5/(г, и, т(гк + 0))

/(г, гк, 9) = /(г, гк, т(гк + 0))+£ ^ ’к ’ К к—(9, - т, (гк + 0)) +

,=1 , 1 1

+2 £ д2/<'а Т + 0)) (91 - т, (к + 0))<е, - т, (г, + 0)),

2 },1=1 691 1

- Ж 5/ (г, 4, т (4 + 0))

/ (г,гк, 9) = / (г,гк,»5(^ + 0)) + £ г кж к—-(9, -т, (гк + 0)) +

,=1 д9,

+'2 £ д2 /^,+ 0)) (9, -т,(г.- + 0))(9, -т,Щ + 0)). (13)

2 , ,1 = 1 , д9,

Рассмотрим числитель в (12). Учитывая (13) и обозначения (11), имеем

I / (г, Ч, 9) р(91Ч + 0)а 9 = / (г, гк, т (гк + 0)) | р(91 р + 0)а 9 +

Ж 5/Мтк + 0))

,=; д9, 0

ь£-

I (91 - т, Ук + 0))р(9 | р + 0)й9 +

I (91 - т, Ук + 0))(9, - т1 Ук + 0))р(9 | гк + 0)й9 :

0

к

к

0

1 Ж 5 / (г,К, т(4 + 0))

= / (г, 4, т (4 + 0))+- £ с,ц + 0). (14)

21,=1 59159,

Аналогично найдем знаменатель в (12):

г - - - 1 Ж 52 /(г, 4, т(4 + 0))

I /(г, 4, 9) р(91 4 + 0)й9 = /(г, 4, т (4 + 0))+- —£ с,, (4 + 0) .(15)

0 2,,,=1 59, 59, ,

Подставляя (14) и (15) в (12), находим

/г (', гк, т(гк + 0)) + 2 £ 5 , + 0)) ^ (гк + 0)

2 , 1=1 59, с9, ----

тг (г) =-----------------------^-----------—1-------------------, г = 1, N. (16)

/о, к, ж,+0))+1 £5 / ^ mе'k+0)) с,, о,+0)

2 ,,1=1 59,се;

Согласно формуле (16), для оценки параметра 9г (9г = тг (г)) в произвольный момент времени г необходимо знать начальные значения т(гк + 0), матрицу кова-

II ||Ж

риаций с(гк + 0) = ||с, (гк + 0)Ц , значения функций / (г, гк,т(гк + 0)) и /(г, гк,т(гк + 0)), г = 1, N . Отметим, что интегрирование при вычислении / (г, гк, т(гк + 0)) и /(г, гк,т(гк + 0)), г = 1, N, производится на конечном отрезке времени [гк, г], г < гк+1.

Начальные значения т (г0 + 0) и С(г0 + 0) задаются согласно апостериорной плотности р(9 | г0 + 0). Если плотность р(9 | г0 + 0) выбрана в соответствии с замечанием к лемме 2, то элементы вектора т (г0 + 0) и элементы матрицы с(г0 + 0) получаются следующим образом:

1 ^ _

»г (г° + 0) = 9(1) - 9(0) I 9гй9г , 1 1 N ,

е<1>

с,, ('0 + 0) = 5(У,,) 9(1) - 9(0) I (9, - т, ('0 + 0))(9, - т, ('0 + 0))й9, , ,,, = ЦУ, (17)

,(0)

где [9г(0), 9(1) ] - отрезок определения равномерной плотности вероятностей

р(9г | г0 + 0) = 1/(9(1) -9(0)) параметра 9г-, г = 1, N, 5(], I) - символ Кронекера.

Рассмотрим формулу (16) для расчета тХ/) в любой момент времени г на полуинтервалах [4, 4+1) (г = 1,N , к = 0, 1, ...). Для вычисления величин тг(г) (г = 1,N), во-первых, необходимо получить выражения для векторов т (4 + 0), матриц ковариаций С(гк + 0) при к = 1, 2, ... . Во-вторых, необходимо получить вторые част-

52 / (г, 4, т (4 + 0)) 52 / (г, 4, т (4 + 0)) —

ные производные функций -----------к------к---, ------г--к----к------, г, ,,I = 1, N .

59,59, 59,59,

В силу (11), аналитические формулы вторых частных производных содержат определенные интегралы, поэтому данные производные рассчитываются численно.

Получим выражения для векторов т (гк + 0) и матриц ковариаций С(гк + 0) при к = 1, 2, ... . Для этого рассмотрим т (гк + 0) (к = 1, 2, ...). С учетом выражений (11) и (8), находим »¿(4 + 0) в виде

19,5(4 - 0,9)р(9 14 - 0)^9

mi (гк + 0) = I 9,р(9 | гк + 0^9 = —-----------------------—,

г к 0 г к I ^к - 0,9)р(9|4 - 0)^9

0

г = , к = 1, 2, ... . (18)

Обозначим (гк - 0,9) = 9is('k - 0,9) и разложим функции s('k - 0,9),

•5 (гк -0,9) (г = 1,N) в ряд Тейлора в окрестности 9 = т(гк -0) с точностью до

о (|л9||2), Д9 = 9 - т (гк - 0). Получим

5(гк - 0, 9) = 5(4 - 0, т(гк - 0)) + £ 5s('k 05т^ 0)) (91 - т, (гк - 0)) +

59,

,=1

N 52

1 N 525(4 - 0, т(4 - 0))

+Т £ к ’ к (9, - т, (4 + 0))(9, - т, (4 + 0)),

- N 55, (4 - 0, т(4 - 0))

5 (гк - 0,9) = 5 (гк - 0 т(гк - 0)) + £-------------59-------(9, - т, (гк - 0)) +

,=1 ^

+2 £ 5 5 ^ Яяе0,Яее('k 0)) (9,- т (гк- 0))(9,- щ (гк- 0)), г =15 N.

2 3,1=1 59 1 ,

Подставляя полученные выражения в (18), получаем

5, ^ - 0т« - »» + 2,£ 152 * ('‘ Яе°,Яе^ - 0)) ст,, (гк + 0) =----------------------- — --------------------------,

/ « ^ 1 ^ 5 5(гк -0,т<гк - 0))

*(4 - 0, т(гк - 0))+- £ ————- с,,

к к 2^ 59,59, 1

],I = ТЖ , к = 1, 2, ..., (19)

где с А = с^к - 0). ______

Получим выражение для с, = Cj/('k - 0) в (18) (],I = 1,N, к = 1, 2, ...). Рассмотрим (11) для ^(1) на полуинтервале времени [4_ь гк) (к = 1, 2, ...). Учитывая фор-

мулу (9) для плотности р(91 г) и обозначения (11), находим

с,, (г) = I (91 - т1 (г))(9, - тг (г)) /(', ^-1,9)р(91^-1 + 0) _ ё9 . (20)

1 0 1 1 I / (г, 4-1,9) рф^к-1 + 0^ 9

0

Обозначим

е, (г, гк-1, 9) = (9, - т, (г))(9, - т1 (г))/(г, 4-1, 9), 4-1 < г < 4, 1, / = ТЖ. (21)

С учетом введенного обозначения, (20) на полуинтервале времени [/к-1, 4) примет вид

| ^ (4 /к _!, 0) р(0|р - + 0)й 0 _

с л (/) = ©-------------------------- , Л, I = 1, N . (22)

| / (/, /к-l, 0) Р(0 1/к-1 + 0)^0

©

Подставляя в (22) / = /к - 0, получаем

| (/к - 0,/к-1,0)Р(0|/к-1 + 0^0 _

л (/к - 0) =©--------------------------------- , ], I = 1, N . (23)

| /(/к- 0, /к-^ 0)Р(01/к-1 + 0)й?0

©

Раскладывая функцию (21) в ряд Тейлора в окрестности точки 0 = т(/к-1 +0) с точностью до о (||Д0||2), Д0 = 0 - т(/к-: + 0), находим

^ (/, /к-1, т (/к-1 + 0))

РА(^/k-l,0) = ра(/,/k-l,т(/к-1 + 0)) + Т-------------¡0----(0“ -ти (/к-1 + 0)) +

и

1 N д2Е]і (1, ^ т(,*-1 + °))

+ О £----------^ *-1-(0и - ти (/*-1 + °))(0, - т, (і*- + °)), у, 1 = 1, N .

2 и~1 д0и 50,

Раскладывая в ряд Тейлора функцию /(/, /к-1,0), стоящую в знаменателе (22), в окрестности точки 0 = т(/к-1 +0) и пренебрегая членами порядка малости

о (||Д0||2), Д0 = 0 - т(/к-1 + 0), получаем

©/ (/,/к-1,0) р(0 |?к-1+0)^0=/(/,/к-1,т (/к-1+0))+2 т9 / ^ тд0/к-1+0)) сл (/к-1+0).

Подставив полученные разложения в (23) и учитывая, что / = /к - 0, получаем (23) в виде

д 2 К,1 ('к- 0 1к-l, т('к-1 + °л

с

1 N д2 Еп (,* - °, I*-1, т (,*-1 + 0))

р,а* -сх,*^,»(<*_1 + °))+- £ ‘ *' ‘1—.

2 и„ 1 д0и

С]1(1* - °) =-

/(/к -0,/к_,т(/м + 0))+2 Тд2/('к-%-двт(-к-1+^

2 иV=1 д0и ^

Л,I = , к = 1, 2, ... , (24)

где сш = с^(4-1 + 0) (и,V = 1, N ).

Формула (24) дает элементы матрицы ковариаций С(/) в момент времени / = /к - 0, необходимые для пересчета апостериорных средних т1 (/к + 0) (i = 1, N) по формуле (19).

Для расчета оценок параметров 0i (/) = т1 (/) (/ = 1, N) на полуинтервале времени [/к, /к+1) по формулам (16) необходимо найти элементы сЛ-(/к + 0) ковариаци-

онной матрицы С(/) в момент времени / = /к + 0 (а,I = 1, N). Имеем

с л (/к + 0) = | (0; - т} (/к + 0))(0/ - т1 (/к + 0)) р(0 | /к + 0)ё 0 , а, I = ЦУ .

©

С учетом введенных обозначений и формулы (8) пересчета плотности вероятностей р(0 | /) в момент времени /к наступления события имеем

I (0; - т} (/к + 0))(0, - т1 (/к + 0))s(/k - 0,0)р(0 | /к - 0)ё0 С л (/к + 0) = © ¡=-------------------------------------------------------------------=-. (25)

15(/к - 0,0)р(0 | /к - 0) ё0

©

Введем обозначение

А а (/к + 0,0) = (0; - т} (/к + 0))(0, - т1 (/к + 0)) 5(/к - 0,0), а, I = .

С учетом введенного обозначения, формула (25) принимает следующий вид:

| Ап (/к + 0,0)р(0|/к - 0)ё0 _

С}1 (/к + 0) = -с----------------------- , А, I = 1, N . (26)

15(/к - 0,0)р(0 | /к - 0)ё0

©

Раскладывая функцию А}1 (/к + 0,0) в ряд Тейлора в окрестности точки

0 = т(/к - 0) с точностью до о (||Д0||2), Д0 = 0 - т(/к - 0), получаем

- * дАп (/к + 0, т(/к - 0))

А А (/к + 0,0) = А А (/к + 0, т (/к - 0)) + Т д0 к----------------(0и - ти (/к - 0)) +

и =1 д0и

1 * д2ап (/к + 0, т(/к - 0)) —

+2 Т-----А кд0 д0 к-(0и - ти (/к - O))(0v - mv (/к - 0)), А, I = 1, N .

2 u,v=1 д0ид0v

Разлагая функцию 5(/к - 0,0) в ряд в точке 0=т(/к - 0), находим

5(/к - 0, 0) = 5(/к - 0, т(/к - 0)) + Т д5(/к 0д0”(/к 0)) (0и - ти (/к - 0)) +

и =1 д0и

+2 Т д 5(/к 0)) (0и- ти (/к- O))(0v- т- (/к- 0)), А, 1 = 1N.

2 u,v =1 д0и д0v

Подставляя данные разложения в (26), получаем формулу для пересчета элементов с?1(/) матрицы ковариации С(/) в момент времени /к наступления события потока:

1 N д2Аа (/к + 0, т(/к - 0))

А}1 (/к + 0, »5(/к - 0)) + - Т--------А кд0 д0 к--------------Сuv (/к - 0)

(. + 0) 2 u,v=1 д0и д0v

СА1 (/к + 0) =

,(,к - 0, т(/к - 0))+2 Т д2^ а,0,д0</к - 0)) Си- (к - 0)

2 и,-=1 д0и д0-

А, 1 = . (7)

4. Алгоритм оценки параметров потока

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

1) в момент времени /0 = 0 исходя из априорных данных о потоке задаются начальные значения оценок параметров потока 0(/0) = т(/0 + 0) и начальная матрица ковариаций С(/0 + 0). В частности, можно считать начальные оценки параметров N независимыми случайными величинами, распределенными равномерно в интервалах допустимых значений параметров потока; в этом случае т(/0 + 0) и

С(/0 + 0) задаются по формулам (17). После задания т(/0 + 0) по формулам (4) вычисляются начальные оценки апостериорных вероятностей ю(ХА- | /0) (А = 1, п), зависящие от т(/0 + 0).

2) Для к = 0 в любой момент времени / (/к < / < /к+1), где /к+1 - момент наблюдения (к + 1)-го события потока, рассчитываются оценки параметров т(/) по формулам (16); вероятности ю(Ха | /) (А = 1,п), входящие в состав формул (16), рассчитываются по формулам (2).

3) В момент времени / = /к+1 наблюдения события асинхронного потока т(/к+1 - 0) рассчитываются по формулам (16) и С(/к+1 - 0) - по формулам (24).

Вероятности ю(Х,- | /) (А = 1,п) - по формулам (2).

4) В момент времени / = /к+1 по формулам (19) и (27) рассчитываются т(/к+1 + 0) и С(/к+1 + 0) соответственно. Вероятности ю(ХА-1 /к+1 + 0) (А = 1,п) пересчитываются по формулам (3).

5) Шаги 2 - 4 алгоритма повторяются для к = 1 и т.д.

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

5. Численный эксперимент

Для численного эксперимента по оценке параметров потока на основе алгоритма, приведенного в разделе 4 настоящей статьи, была реализована программа расчета оценок на языке С++.

Первый этап расчета - имитационное моделирование потока с заданными параметрами.

Второй этап расчета - вычисление оценок параметров, а также апостериорных вероятностей и элементов матрицы ковариации С по формулам (2) - (4), (16), (17), (19), (24), (27).

В приведенном ниже примере расчеты были произведены для следующих значений параметров: п =2 (асинхронный поток имеет два состояния), Х1 = 0,5, Х2 = 0,05, а1 = а12 = 0,08, а2 = а21 = 0,04. При этом начальная плотность распределения параметров потока выбрана равномерной: параметр Х1 равномерно распре-

делен в интервале [0,337;0,836], параметр Х2 - в интервале [0,021;0,073], параметр а! - в интервале [0,04;0,091], параметр а2 - в интервале [0,031;0,048]. Время моделирования Т изменяется от 0 до 100 ед. времени.

На рис. 1 - 4 приведены траектории изменения оценок параметров потока в зависимости от времени моделирования Т.

Рис. 1. Траектория оценки параметра Х1(г)

^2(0

0,052

0,050

0 20 40 60 80 Т

Рис. 2. Траектория оценки параметра Х2(г)

«1(01

0,08---------------------------------------------------------------

0,07

0,06

20 40 60

Рис. 3. Траектория оценки параметра а1(/)

Т

СС2О)

0,0400

0,0398

0 20 40 60 80 Т

Рис. 4. Траектория оценки параметра а2(г)

Скачки траекторий соответствуют моментам наступления событий потока - в эти моменты оценки параметров пересчитываются по формулам (19) и претерпевают разрывы первого рода.

Анализ приведенного в настоящей статье, а также многочисленных проведенных численных экспериментов показывает, что оценки параметров потока имеют достаточно хорошую стабильность, в том числе в случае равномерного распределения начальных оценок параметров в интервалах, допускающих отклонение от истинного значения параметра на 100 %. Колебание оценок для

Я,1(/) происходит в диапазоне от 0,578281 до 0,737368; для Х2(?) - в диапазоне

от 0,04702 до 0,0486; для ос1(/) - в диапазоне от 0,06279 до 0,06572; для а2(/) -в диапазоне от 0,0396 до 0,03968. Еще раз подчеркнем, что полученные оценки обеспечивают минимум среднеквадратического отклонения оценок параметров от истинных значений.

Заключение

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

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

ЛИТЕРАТУРА

1. Башарин Г.П., Кокотушкин В.А., Наумов В.А. О методе эквивалентных замен расчета

фрагментов сетей связи. Ч.1 // Изв. АН СССР. Техн. кибернетика. 1979. № 6. С. 92-99.

2. Башарин Г.П., Кокотушкин В.А., Наумов В.А. О методе эквивалентных замен расчета фрагментов сетей связи. Ч.2 // Изв. АН СССР. Техн. кибернетика. 1980. № 1. С. 55-61.

3. NeutsM.F. A versatile Markov point process // J. Appl. Probab. 1979. V. 16. P. 764-779.

4. Горцев А.М., Нежельская Л.А. Оценка параметров синхронного альтернирующего пуас-соновского потока событий методом моментов // Радиотехника. 1995. № 7-8. С. 6-10.

5. Горцев А.М., Шмырин И.С. Оптимальная оценка параметров дважды стохастического пуассоновского потока событий при наличии ошибок в измерениях моментов наступления событий // Изв. вузов. Физика. 1999. № 4. С. 19-27.

6. Бушланов И.В., Горцев А.М., Нежельская Л.А. Оценка параметров синхронного дважды стохастического потока событий // Автоматика и телемеханика. 2008. № 9. С. 76-93.

7. Горцев А.М., Зуевич В.Л. Оптимальная оценка состояний асинхронного дважды стохастического потока событий с произвольным числом состояний // Вестник Томского государственного университета. Управление, вычислительная техника и информатика. 2010. № 2(11). С. 44-65.

8. Хазен Э.М. Методы оптимальных статистических решений и задачи оптимального управления. М.: Сов. радио, 1968. 256 с.

Горцев Александр Михайлович

Зуевич Владимир Леонидович

Томский государственный университет

E-mail: [email protected]; [email protected] Поступила в редакцию 23 августа 2011 г.