Имитация марковского процесса

4.6.1. Моделирование дискретной цепи Маркова. Рассмотрим дискретную цепь Маркова или марковский процесс с дискретным временем перехода из одного состояния в другое. Математическая модель дискретной цепи Маркова задается матрицей вероятностей переходов ||Р_ij||, имеющей вид

, (4.12)

где Р_ij - вероятность перехода из состояния z_i в состояние z_j в некоторый дискретный момент времени kDt (k=1,2,3,...).

Начальное состояние марковского процесса определяется матрицей-строкой начальных вероятностей ||Р₀||=|Р₁(0), Р₂(0), ..., Р_n(0)|,гдеР_i(0) - вероятность нахождения процесса в z_i-м состоянии при t=0. Вероятности перехода Р_ij не зависят от времени, т.е. процесс является однородным.

На рис. 4.23 приведен обобщенный алгоритм имитации дискретной цепи Маркова.

 

 

Рис. 4.23

 

Подпрограммы WWOD и WIWOD реализуют интерфейсную часть имитационной программы. В подпрограмме WWOD осуществляется ввод элементов массива Р0[I], в который заносятся вероятности Р_i(0), и массива Р[I,J], в который заносятся вероятности Р_ij. Определяется начальный такт моделирования T=0 и заданное число тактов моделирования TZ.

Подпрограмма OРRZ0 предназначена для определения начального состояния, а подпрограмма OРRZ – для определения состояний в процессе моделирования смены состояний. Подпрограмма STAT предназначена для набора статистических данных.

Моделирование марковского процесса основано на принципе имитации системы случайных событий.

Определяем числовые границы

.

Исходя из значения случайного числа Р₀, генерированного датчиком случайных чисел, определяется номер rначального состояния z_r(0), для которого будет справедливо условие

.

Затем датчик случайных чисел вырабатывает случайное число Р₁, которое также сравнивается с границами

.

Путем сравнения устанавливается очередное состояние и подобным образом осуществляется моделирование дальше.

Рассмотрим реализации подпрограмм.

На рис. 4.24 приведен алгоритм подпрограммы OРRZ0.

На рис. 4.25 приведен алгоритм подпрограммы OРRZ.

На рис. 4.26 приведен алгоритм подпрограммы STAT.

В алгоритме подпрограммы OРRZ0 (см. рис. 4.24) в блоке 1 вырабатывается число Р датчиком случайных чисел.

Затем реализуется цикл по переменной I для сравнения числа Р с элементами массива Р0[I]. Для этого введен идентификатор А, который в блоке 2 определен А=0.

 

Рис. 4.24

 

 

Рис. 4.25

 

 

Рис. 4.26

 

При наращивании переменной I идентификатор А будет принимать последовательно значения: I=1, А=Р₁(0); I=2, А=Р₁(0)+Р₂(0); I=3, А=Р₁(0)+Р₂(0)+Р₃(0); …; I=n, А=Р₁(0)+Р₂(0)+Р₃(0)+…+Р_n(0).

При первом же выполнении условия Р£A считается, что найден индекс начального состояния z_i(0). Таким образом, выходным параметром подпрограммы OРRZ0 является значение индекса I, при котором выполнено условие Р£А=Р₁(0)+Р₂(0)+Р₃(0)+…+Р_i(0).

В алгоритме подпрограммы OРRZ (см. рис. 4.25) в блоке 1 датчиком случайных чисел формируется случайное число РÎ[0,1]. Затем реализуется цикл по переменной J для сравнения числа Р с элементами массива Р[I,J]. Введен идентификатор В. При наращивании переменной J идентификатор В будет принимать последовательно значения: J=1, А=Р(I,1); J=2, А=Р(I,1)+Р(I,2); …; J=n, А=Р(I,1)+Р(I,2)+Р(I,3)+…+Р(I,n).

При выполнении условия Р£B считается, что найдет индекс состояния z_j(T) в текущем такте моделирования T.

Выходным параметром подпрограммы OРRZ является значение индекса J, при котором выполнено условие Р£В=Р(I,1)+Р(I,2)+Р(I,3)+…+Р(I,J).

В алгоритме подпрограммы STAT (см. рис. 4.26) в блоке 1 в счетчиках K[J] осуществляется подсчет частот появления событий z_j. Затем в блоке 2 осуществляется присвоение значения индекса J индексу I, т.е. на следующем такте моделирования в подпрограмме OРRZ будет выполнен анализ j-й сроки матрицы вероятностей переходов||Р_ij||. Алгоритм имитации дискретной цепи Маркова может быть представлен в более сокращенном виде, как показано на рис. 4.27. Объединим матрицы ||Р₀|| и ||Р_ij|| в обобщенную матрицу вида

. (4.13)

В подпрограмме WWOD осуществляется ввод элементов массива Р[I,J], в который заносятся вероятности Р_ij, причем , . Определяется начальное значение индекса I=0, начальный такт моделирования T=0 и заданное число тактов моделирования TZ.

Так как в первом такте моделирования (Т=1) I=0, то в первом такте будет рассматриваться верхняя строка обобщенной матрицы ||Р_ij||, т.е. фактически матрицы ||Р₀||.

В блоке 3 формируется случайное число РÎ[0,1]. Затем реализуется цикл по переменной J для сравнения числа Р с элементами массива Р[I,J] при I=0 (см. блоки 4 – 6). При выполнении условия Р£В в соответствующий счетчик K[J] добавляется единица (см. блок 7). В блоке 8 индексу I присваиваются значения индекса J. Процесс моделирования продолжается до тех пор, пока не будет выполнено условие окончания моделированияT<TZ.

 

 

Рис. 4.27

 

В счетчиках K[J] после окончания моделирования будут подсчитаны частоты пребывания марковского процесса в состояниях z_j, . В подпрограмме WIWOD осуществляется вывод значений счетчиков K[J], а также вывод эмпирических оценок финальных вероятностей марковского процесса, определяемых по формуле

.

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

- распределением вероятностей начального состояния процесса Р_i(0) в момент t₀;

- матрицей вероятностей переходов ||Р_ij|| (12);

- матрицей-строкой функций распределения времени пребывания в состоянии |A|={A₁(t), A₂(t), …, A_n(t)}, где A_i(t)- функция распределения времени пребывания марковского процесса в состоянии z_i.

На рис. 4.28 приведен обобщенный алгоритм имитации вложенной цепи Маркова.

 

 

Рис. 4.28

 

Подпрограммы WWOD и WIWOD предназначены для реализации интерфейса имитационной программы. В подпрограмме WWOD осуществляется ввод элементов массива Р0[I], массива Р[I,J], а также параметров функций распределения A_i(t). Определяется T=0 и TZ. В подпрограмме OРRZ0 (см. блок 3) определяется начальное состояние, в подпрограмме OРRZ определяются состояния в процессе моделирования смены состояний.

Подпрограмма STAT1 предназначена для набора статистических данных о частотах пребывания марковского процесса в состояниях.

В подпрограмме OРRTAU определяется время пребывания марковского процесса в состоянии, определенном в подпрограмме OРRZ. Подпрограмма STAT2 предназначена для набора статистических данных о времени пребывания марковского процесса в состояниях.

Рассмотрим реализацию подпрограмм OРRZ0, OРRZ, OРRTAU.

Подпрограммы OРRZ0, OРRZ имеют такой же вид, как и аналогичные подпрограммы для дискретной цепи Маркова (см. рис. 4.24 и рис. 4.25). Выходным параметром подпрограммы OРRZ является индекс J состояния z_j(T) в текущем такте моделирования T.

Схема алгоритма подпрограммы OРRTAU будет иметь вид алгоритмов генерации случайной величины, рассмотренных в разд. 4.5, метод обратных функций, метод ступенчатой аппроксимации, использование предельных теорем. О схемной реализации подпрограмм STAT1 и STAT2 будет сказано ниже.

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

- распределением вероятностей состояния Р_i(0) в момент t₀;

- матрице L, имеющей вид

, (4.14)

где l_ij -интенсивностей переходов марковского процесса из состояния z_i в состояние z_j.

Способ моделирования сводится к следующему.

Пусть в момент t состояние марковского процесса z(t)=z_i.Известен момент времениt_i, в который марковский процесс должен покинуть состояние z_i. Для определения следующего состояния z_i+1 и момента t_i+1 перехода в последующее состояние генерируется ряд независимых случайных величин t_ijпо показательному закону распределения с параметрами i-й строки матрицы ||l_ij ||. Затем определяется минимальное значение из полученной совокупности независимых случайных величин t_ij, т.е.

. (4.15)

Следующий момент изменения состояния будет t_i+1=t_i+t_ij, а состояние, в которое переходит система в момент t_i+1, будет z_j, где z_j - состояние, при котором t_zjбыло минимальным из всех значений совокупности независимых случайных величин {t_ii}. Физическая интерпретация такого подхода к моделированию марковского процесса с непрерывным временем позволяет разрабатывать имитационные модели реальных систем.

Схема алгоритма моделирования вложенной цепи Маркова по данному способу приведена на рис. 4.29.

 

 

Рис. 4.29

 

Отличие алгоритма по схеме развернутой рекуррентной имитации (см. рис. 4.29) от алгоритма, приведенного на рис. 4.28, состоит в блоках 4 и 5, в которых реализована генерация случайных величин t_kiи формулы (4.15).

На рис. 4.30 приведена схема алгоритма подпрограммы OРRTAU для данного способа. На рис. 4.31 приведена схема алгоритма подпрограммы MINTAU.

 

 

Рис. 4.30

 

В подпрограмме OРRTAU (см. рис. 4.30) входным параметром является индекс I предыдущего состояния марковского процесса.

Для каждого значения J в блоке 4 по методу обратных функций определяется значение случайной величины TAU[J], отождествляющей собой описанную выше величину t_ii. Таким образом, в данной подпрограмме будет сформирован массив TAU[J].

В подпрограмме MINTAU (см. рис. 31) входным параметром являются элементы массива TAU[J], из которых необходимо выбрать минимальный элемент.

В блоках 1 – 5 осуществляется выбор минимальной величины (MIN) из массива TAU[J]. Происходит это следующим образом. При переменной массива J=1 первый элемент массива TAU[1] определяется, как минимальный элемент MIN (см. блок 1).

 

 

Рис. 4.31

 

Затем второй элемент массива TAU[2] сравнивается со значением MIN (см. блок 3). Если второй элемент массива TAU[2] меньше значения MIN, то идентификатору MIN присваивается значение TAU[2] (см. блок 4). Если условие TAU[2]<MIN не выполняется, то рассматривается условие TAU[3]<MIN для элемента массива TAU[3] и т.д.

В блоках 6 – 8 определяется индекс J элемента массива TAU[J], равного числу MIN. Это индекс J определяет индекс состояния z_J, в который переходит марковский процесс.