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

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

Для начала, вычислим матрицу дисперсии D, по формуле D=N*(2*Ndg – E)-Nsq, для этого нам понадобятся матрицы:

Ndg – матрица N, в которой все элементы помимо главной диагонали =0:

1,607652	0	0	0
0	1,343101	0	0
0	0	1,125356	0
0	0	0	1,282051

 

Nsq – матрица N, в которой каждый элемент возведен в квадрат:

2,584544	1,035308	0,392854	0,507301
0,014908	1,803921	0,182628	0,883921
0,103382	0,041412	1,266426	0,020292
0,001342	0,162353	0,016437	1,643655

 

E – Единичная матрица.

В итоге, матрица D равна:

0,976892	0,680405	0,391068	0,606732
0,25558	0,46082	0,351864	0,586602
0,608905	0,30173	0,14107	0,202515
0,079805	0,517069	0,143911	0,361604

 

 

Теперь, когда мы имеем матрицу дисперсии D, нужно найти среднеквадратичное отклонение числа пребывания процесса от среднего G, для этого нам нужно вычислить матрицу, каждый элемент которой получен путем извлечения квадратного корня из соответствующего элемента матрицы D:

0,988379	0,824866	0,625355	0,77893
0,50555	0,678837	0,593181	0,7659
0,780324	0,5493	0,375593	0,450016
0,282498	0,719075	0,379356	0,601335

 

Теперь найдем среднеквадратичное отклонения трудоемкости процесса от среднего, так: Gср=G11*C1+ G12*C2+ G13*C3+ G14*C4=14,00535.

 

                      Оценка предельных вероятностей пребывания процесса во множестве невозвратных состояний от среднего:

Для оценки предельных вероятностей пребывания процесса в множестве T~, воспользуемся 2-мя методами:

1. Возведем матрицу P в предельно большую степень. limPk при k-> к бесконечности, и получим:

0	0	0	0	0,286325	0,286325	0,42735
0	0	0	0	0,217949	0,217949	0,564103
0	0	0	0	0,457265	0,457265	0,08547
0	0	0	0	0,115385	0,115385	0,769231
0	0	0	0	0,5	0,5	0
0	0	0	0	0,5	0,5	0
0	0	0	0	0	0	1

 

Как видно, матрица Q стремится к нулевой матрице.

2. Для проверки правильности результата, воспользуемся спектральным разложением матрицы P. Когда матрица P возводится в k степень, она имеет вид

Qk	R(k)
0	Wk

При предельном возведении, матрица Q->0, а матрицы R(k) и Wk, принимают вид:

R(k):

0,286325	0,286325	0,42735
0,217949	0,217949	0,564103
0,457265	0,457265	0,08547
0,115385	0,115385	0,769231

 

 

Wk:

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

 

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

0

0

0

0

0,217949

0,217949

0,564103

 

5. Список литературы:

• Г.А Доррер «Методы моделирования дискретных систем»   Красноярск 2004.

• http://ru.wikipedia.org – статья «Введение в Цепи Маркова».
• http://www.statsoft.ru – статья «Применение цепей Маркова».

6. Приложения:

 

//---------------------------------------------------------------------------

 

#include <vcl.h>

#include <stdlib.h>

#include <stdio.h>

#pragma hdrstop

#define siz 7

#define siz1 425

#define siz2 425

 

#include "Unit1.h"

//---------------------------------------------------------------------------

#pragma package(smart_init)

#pragma resource "*.dfm"

TForm1 *Form1;

int matr[siz][siz]; //матрица.

double met[25][siz]; //массив графика.

double srmet[siz];   //Массив средних значений.

AnsiString strk[25]; //Строки статистики.

int vekt[siz]; //Вектор, выбора пути.

int vind[siz]; //Индексы соответствующих элементов вектора vekt.

int h=0;

int pos1=0; //Переменная позиции, запоминает позицию выбранного.

int x=1000; //Количество тестов.

int y=15;

int dood=0;

//---------------------------------------------------------------------------

__fastcall TForm1::TForm1(TComponent* Owner)

        : TForm(Owner)

{

}

//---------------------------------------------------------------------------

void __fastcall TForm1::FormShow(TObject *Sender)

{

int i,j;

 

//РИСУЕМ  МЕТКИ//

//Метка

AnsiString sf="";

for(i=0;i<440;i=i+40)

{

if(i<400) sf=i/40;

else sf=i/400;

img1->Canvas->Pen->Color=(TColor)RGB(0,140,140);//Рисуем рамку

img1->Canvas->Pen->Width=1;

if(i<400) img1->Canvas->TextOutA(2,siz2-5-i-12,"0."+sf);

else  img1->Canvas->TextOutA(2,siz2-5-i-12,sf);

 

}

 

//-----

//************//

 

//МОЯ МАТРИЦА//

matr[0][0]=3;matr[0][1]=5;matr[0][2]=2;

matr[1][2]=3;matr[1][3]=7;

matr[2][0]=2;matr[2][5]=8;

matr[3][1]=3;matr[3][4]=1;matr[3][6]=6;

matr[4][4]=5;matr[4][5]=5;

matr[5][4]=5;matr[5][5]=5;

matr[6][6]=10;

//===========//

for(i=1;i<siz+1;i++) //Заполнение матрицы элементами массива.

{

        for(j=1;j<siz+1;j++)

        {

        matr1->Cells[i][j]=matr[j-1][i-1];

        }

}

for(i=1;i<siz+1;i++)//Нумеруем строки.

{

matr1->Cells[0][i]=i-1;

matr1->Cells[i][0]=i-1;

}

for(i=1;i<siz+1;i++)//Нумеруем строки.

{

matr2->Cells[i][0]=i-1;

}

for(i=1;i<16;i++)//Нумеруем строки.

{

matr2->Cells[0][i]=i;

}

for(i=0;i<7;i++)

{

if(i==0) matr4->Cells[i][0]=1;

else  matr4->Cells[i][0]=0;

}

 

}

//---------------------------------------------------------------------------

void __fastcall TForm1::start1Click(TObject *Sender)

{

int i,j,g;

if(test1->Text!="") x=test1->Text.ToInt();

else x=100;

y=shag1->Text.ToInt();

srand(time(NULL));

AnsiString st="";

//=============Обновление  массива и матрицы==============//

for(i=0;i<siz;i++)//Таблицу в массив.

{

        for(j=0;j<siz;j++)

        {

        st=matr1->Cells[j+1][i+1];

        if(st!="") matr[i][j]=st.ToDouble();

        }

}

for(i=0;i<siz;i++) //Массив в таблицу.

{

        for(j=0;j<siz;j++)

        {

        matr1->Cells[j+1][i+1]=matr[i][j];

        }

}

//====================Начало  тестов==========================//

for(g=0;g<x;g++) //***

{

for(i=0;i<siz;i++)//Заполняем начальный вектор.

{

vekt[i]=matr[1][i];

vind[i]=i;

}

//=====================15 шагов==========================//

for(h=0;h<y;h++) //**

{

//Упорядочим  вектор.

int k=100,l=-1;

for(i=0;i<siz;i++)

{

   for(j=i;j<siz;j++)

   {

        if(vekt[j]<k && vekt[j]!=0)

        {

        k=vekt[j];

        l=j;

        }

   }

   if(l>-1)

   {

        k=vekt[i];  //Ставим минимальный на 1 место.

        vekt[i]=vekt[l];

        vekt[l]=k;

        k=vind[i];

        vind[i]=vind[l]; //И ставим его индекс.

        vind[l]=k;

   }

   k=11;

   l=-1;

}

//Генерация  числа.

k=rand()%100;   //Генерируем число от 0 до 9.

k=k%10;

i=0;

int d=0;

while(k>(vekt[i]-1+d))  //Выбираем по какому пути идем.

{

        d=d+vekt[i];

        i++;

}

pos1=vind[i]; //Фиксируем индекс выбранного.

met[h][pos1]++;  //Увеличиваем статистическую информацию на шаге по вершине.

for(i=0;i<siz;i++)   //Записываем в вектор новую строку.

{

vekt[i]=matr[pos1][i];

vind[i]=i;

}

}

 

}

double sk=0;

AnsiString re1="";

//**********Считаем среднее**********//

for(i=0;i<siz;i++)

{

        for(j=0;j<y;j++)

        {

             sk=sk+met[j][i];

        }

        re1=sk/x;

        re1.Delete(6,re1.Length());

        srmet[i]=re1.ToDouble();

        teh1->Items->Add(srmet[i]);

        sk=0;

        if(i==1) srmet[i]=srmet[i]+1;

        matr3->Cells[i][0]=srmet[i];

}

//***********************************//

re1="";

for(i=0;i<y;i++)//Расчет как в икселе.

{

        for(j=0;j<siz;j++)//Высчитываем новые вероятности.

        {

        met[i][j]=met[i][j]/x;

        re1=met[i][j];

        re1.Delete(6,re1.Length());

        met[i][j]=re1.ToDouble();

        }

}

for(i=0;i<y;i++)//вывод статистической информации.

{

        for(j=0;j<siz;j++)

        {

                matr2->Cells[j+1][i+1]=met[i][j];

        }

}

}

//---------------------------------------------------------------------------

//Очистка статистики.

void __fastcall TForm1::ClearClick(TObject *Sender)

{

int i,j;

for(i=0;i<y;i++)//Очистка статистическоф информации.

{

        for(j=0;j<siz;j++)

        {

                met[i][j]=0;

                matr2->Cells[j+1][i+1]="";

        }

}

}

//---------------------------------------------------------------------------

 

 

 

void __fastcall TForm1::gr1Click(TObject *Sender)

{

int j=0,a=0;

double i=0,k=20,l=siz2;

//**//

if(dood==0) img1->Picture=NULL;

img1->Canvas->Pen->Color=(TColor)RGB(23,23,150);

img1->Canvas->Pen->Width=4;

img1->Canvas->MoveTo(0,1);

img1->Canvas->LineTo(0,siz2);

img1->Canvas->MoveTo(0,siz2);

img1->Canvas->LineTo(siz1,siz2);

img1->Canvas->MoveTo(siz1,siz2);

img1->Canvas->LineTo(siz1,0);

img1->Canvas->MoveTo(siz1,0);

img1->Canvas->LineTo(0,0);

img1->Canvas->Pen->Color=(TColor)RGB(0,0,0);

AnsiString sf="";

for(i=0;i<440;i=i+40)

{

if(i<400) sf=i/40;

else sf=i/400;

if(i<400) img1->Canvas->TextOutA(2,siz2-5-i-12,"0."+sf);

else  img1->Canvas->TextOutA(2,siz2-5-i-12,sf);

}

 

//-----

//**//

img1->Canvas->Pen->Color=(TColor)RGB(250,0,0);

img1->Canvas->Pen->Width=2;

img1->Canvas->MoveTo(k,siz2-5);

for(j=0;j<siz;j++)//По ячейкам строки.

{

i=met[0][j]*100;

l=siz2-5-i*4;

if(j!=0) img1->Canvas->LineTo(k,l);

img1->Canvas->MoveTo(k,l);

k=k+40;

}

}