Моделирование динамики систем на основе цепей Маркова с дискретным временем

 

 

 

 

 

Pk=(UɅU)k=UɅkU-1

Причем Ʌ k

ʎk 1	….	0
….	….	….
0	….	ʎk n

 

Таким образом, динамику изменения матрицы Рк легко оценить по поведению ʎi,   i = 1,...,п .

Мы уже упоминали, что матрица Р, определяющая однородную цепь Маркова, является стохастической матрицей. Особенностью этой матрицы является то, что ее максимальное собственное число равно 1, и ему соответствует собственный вектор, составленный из единиц: [l,l,...l].

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

Xпред= limX(tk) k→∞.

Напомним, что вектор состояния X{tk) определяете» формулой:

X(tk)=X(t0)Pk

а к-я степень канонической стохастической матрицы имеет структуру: Pk=

Qk	R(k)
0	Wk

 

При большом числе шагов, как мы установили, матрица Qk  -» 0, а матрицы Wk и R(k) стремятся к некоторым предельным Wk -» Wnред, R(k) -» Rпред. Таким образом, в предельном случае: Рпред=limPk , k→∞, =

 

0	Rпред
0	Wпред

 

Матрицу Рпред можно найти описанным выше способом,

выполнив спектральное разложение матрицы Р и положив k→∞. Тогда

Xпред= X(t0)*Pпред

Однако существует более простой способ определении Хпред. Он основан на том, что при большом числе шагов вектор X(tk) сходится к Хпред . Из соотношения

X(tk+1) = X{tk)*P, при k→∞,, следует, что

X(tk) -»X{tk+])-»Xпред

и

Xпред=Xпред*P

Решая уравнение относительно неизвестных x1пред, x2пред,   …, xnпред   при дополнительном условии

Σxi пред=1, при i→∞,. Получим вектор Xпред.

3. Постановка задачи и исходные данные: Задан граф цепи Маркова с начальным состоянием 2:

 

 

 

 

 

 

 

 

Матрицп переходных вероятностей для данного графа показана ниже:

	1	2	3	4	5	6	7
1	0,3	0,5	0,2	0	0	0	0
2	0	0	0,3	0,7	0	0	0
3	0,2	0	0	0	0	0,8	0
4	0	0,3	0	0	0,1	0	0,6
5	0	0	0	0	0,5	0,5	0
6	0	0	0	0	0,5	0,5	0
7	0	0	0	0	0	0	1

 

 

Начальный вектор Xt0:

0

1

0

0

0

0

0

 

4. Практическая часть:

            Переходные вектора и график полученные теоретически и практически:

Для 15 шагов были высчитаны вектора переходных вероятностей Xtk , а так же были построены графики к ним, все вычисления велись в Microsoft Excel:

Шаг 1: Xt1

0

0

0,3

0,7

0

0

0

      

 

 

 

 

 

 

Шаг 2: Xt2

0,06

0,21

0

0

0,07

0,24

0,42

      

 

Шаг 3: Xt3

0,018

0,03

0,075

0,147

0,155

0,155

0,42

      

Шаг 4: Xt4

0,0204

0,0531

0,0126

0,021

0,1697

0,215

0,5082

      

 

 

Шаг 5: Xt5

0,00864

0,0165

0,02001

0,03717

0,19445

0,20243

0,5208

       

Шаг 6: Xt6

0,006594

0,015471

0,006678

0,01155

0,202157

0,214448

0,543102

        

Шаг 7: Xt7

0,0033138

0,006762

0,0059601

0,01083

0,209458

0,213645

0,550032

      

 

 

 

 

Шаг 8: Xt8

0,00218616

0,004906

0,00269136

0,004733

0,212634

0,216319

0,55653

      

Шаг 9: Xt9

0,00119412

0,002513

0,001908975

0,003434

0,21495

0,21663

0,55937

      

Шаг 10: Xt10

0,000740031

0,001627

0,000992754

0,001759

0,216133

0,217317

0,56143

      

 

 

Шаг 11: Xt11

0,00042056

0,000898

0,00063619

0,001139

0,216901

0,217519

0,562486

      

Шаг 12: Xt12

0,000253406

0,000552

0,000353442

0,000628

0,217324

0,217719

0,563169

      

Шаг 13: Xt13

0,00014671

0,000315

0,000216284

0,000386

0,217585

0,217804

0,563546

      

 

Шаг 14: Xt14

8,72698E-05

0,000189

0,000123912

0,000221

0,217733

0,217868

0,563778

      

Шаг 15: Xt15

5,09634E-05

0,00011

7,42371E-05

0,000132

0,217822

0,217899

0,563911

      

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Далее представлена таблица векторов, полученная программным путем:

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

Шаг 1:

Шаг 2:

Шаг 3:

 

 

 

 

 

 

 

Шаг 4:

Шаг 5:

 

 

 

 

 

 

 

Шаг 6:

Шаг 7:

 

 

 

 

 

 

 

Шаг 8:

Шаг 9:

 

 

 

 

 

 

 

Шаг 10:

Шаг 11:

 

 

 

 

 

 

 

Шаг 12:

Шаг 13:

 

 

 

 

 

 

 

Шаг 14:

Шаг 15:

 

            

 

 

 

 

 

     Структурирование матрицы P, и выделение подмножеств T и T~:

Структурируем матрицу P, и выделяем из нее матрицы Q,R,W, данные матрицы представлены ниже:

Q:

0,3	0,5	0,2	0
0	0	0,3	0,7
0,2	0	0	0
0	0,3	0	0

 

R:

0	0	0
0	0	0
0	0,8	0
0,1	0	0,6

 

W:

0,5	0,5	0
0,5	0,5	0
0	0	1

После этого нужно выделить множество невозвратных состояний, оно включает в себя T={S1,S2,S3,S4}, и эргодическое множество которое состоит изT~={S5,S6,S7}. Отметим что состояние S7 поглощающее, то есть попав в него процесс его уже не покинет.

 

    Сравнение теоретически полученной матрицы N с результатом работы программы.                 Вычисление средней трудоемкости Cср:

Теперь необходимо вычислить матрицу N, которая показывает среднее число тактов пребывания процесса в каждом из невозвратных состояний. Матрица N вычисляется по формуле N=(E-Q)-1, где E – единичная матрица.

Q:	0,3	0,5	0,2	0
	0	0	0,3	0,7
	0,2	0	0	0
	0	0,3	0	0
E:	1	0	0	0
	0	1	0	0
	0	0	1	0
	0	0	0	1

 

 

 

 

 

Матрица N:

1,607652	1,017501	0,626781	0,712251
0,1221	1,343101	0,42735	0,940171
0,32153	0,2035	1,125356	0,14245
0,03663	0,40293	0,128205	1,282051

Каждая строка матрицы N показывает среднее число тактов пребывания процесса во множестве невозвратных состояний. В соответствии с вариантом заданы трудоемкости:

3

3

6

9

В нашем примере, старт произошел из состояния S2, следовательно помножим каждую трудоемкость на соответствующий элемент строки 2 матрицы N, и получим:

0,3663

4,029304

2,564103

8,461538

В сумме это Cср=15,42125.

Если сравнить 2 строку матрицы N, с полученным программно результатом, то видно что они практически идентичны: