西安电子科技大学(2016年度)随机过程与排队论班级:XXXXXXX姓名:XXX XXX学号:XXXXXXXXXX XXXXXXXXXXX一步转移概率矩阵收敛快慢的影响因素作者姓名:XXX XXX 指导老师姓名:XXX(西安电子科技大学计算机学院,陕西西安)摘要:根据课程教材《排队现象的建模、解析与模拟【西安电子科技大学出版社曾勇版】》,第[1.3马尔可夫过程]中,马尔可夫过程链n时刻的k步转移概率结果,当k=1时,得到一步转移概率。
进而得到一步转移概率矩阵P(1)。
为研究此一步转移概率矩阵(下称一步矩阵)的收敛特性以及影响其收敛快慢的因素,使用MATLAB实验工具进行仿真,先从特殊矩阵开始做起,发现规律,然后向普通矩阵进行拓展猜想,并根据算术理论分析进行论证,最终得出一步矩阵收敛快慢的影响因素。
关键词:一步转移概率矩阵 MATLAB 仿真猜想一、问题概述我们讨论时一步矩阵的特性应从以下两方面来分析:(1)矩阵P(n)在满足什么条件时具有收敛特性;对于矩阵P(n),当P(n)=P(n+1)时,我们说此矩阵具有收敛特性,简称矩阵P(n)收敛。
(2)若一个一步矩阵具有收敛特性,那么其收敛速度与什么有关?首先,我们需要明确什么是一步矩阵收敛:对于一般的一步矩阵P 、矩阵An+1、矩阵An,若有: An+1=AnP=An那么称该一步转移矩阵可收敛。
二、仿真实验1、仿真环境本次采用的是MATLAB仿真实验软件进行仿真实验2、结果与分析【1】、特殊矩阵:单位矩阵与类单位矩阵从图(1)和图(2)可以看出,单位矩阵不具有收敛特性,类单位矩阵并非单位矩阵但是经过n次后也变为单位矩阵,所以此矩阵也不具有收敛特性。
此类矩阵也易证明其不具有收敛性。
图(1)单位矩阵图(2):类单位矩阵【2】、一般单位矩阵图(3):一般一步矩阵Ⅰ图(4):一般一步矩阵从图(3)和()可以看出他们分别在18次和4次后收敛到一个稳定的值3、根据实验的猜想根据在单位矩阵和一般单位矩阵和一般一步矩阵中得到的结果,可以对得出如下结论:类单位矩阵、单位矩阵是不具有收敛性的,而一般的一步矩阵是有收敛性的,而且收敛速率有快有慢。
对于上面结论中的状况,我们首先观察如上四个矩阵,不难发现,在矩阵收敛的最终结果矩阵中,其每行和均为1,而且每列上的值均为相同值。
最终概率分布结果也是矩阵收敛后的一行。
所以根据上述的结果及分析做出如下猜想:每一列比较均匀的矩阵收敛速度较快;与类单位矩阵类似的矩阵收敛速度较慢。
在极限情况下,有如下情况:[1]、列相同矩阵已经是收敛矩阵[2]、已经是是类单位矩阵的,不会收敛。
下面是刻画矩阵收敛速度的方法:根据矩阵的行列式的值,当矩阵的行列式的绝对值为1时,矩阵为类单位矩阵,不会收敛,是收敛最慢的极限。
当矩阵行列式为0时,是收敛最快的极限。
根据以上分析,行列式值越接近1,越与类单位矩阵类似,稳定速率越慢。
矩阵的行列式值越接近0,收敛越快。
作为例证,我们计算一下上述两个矩阵的行列式的值:图(5):矩阵B的行列式值图(6):矩阵C的行列式值从上述的验证中可以看到矩阵B的行列式的绝对值为0.0255 而矩阵C的行列式绝对值为6*10^(-6),远小于行列式B中的值,而正好矩阵B的列值相似度要小于矩阵C。
4、分析与证明我们先看类单位矩阵的行列式的值 为1 而且不难证明所以得一步转移概率矩阵的行列式的值得绝对值都在[0,1]之间。
假设一个n 阶一步转移概率矩阵Q 其行列式的表达式为:Det(Q)=a11*(-1)1+1Det(c (11))+a12*(-1)1+2Det(c (12))….+a1n*(-1)1+nDet(c(1n))由上式可以看出,若列值的差值越大,那么行列式的值就取决于该列的值中的较大的值,若行列式的列差值比较小,那么最终行列式降阶到2阶时,计算得到的值为对角线相减,由于列值相差小,所以所得到的值也会相对较小,也会比较靠近0。
而差值越大,决定因素也会由列中较大值决定以此类推,到最后降阶到2阶时起决定因素的系数都为列中的较大值,而最后的二阶行列式由于差值较大所以计算的结果也会比较大,整体行列式的值都会靠近1。
换个角度 可以将单位矩阵看成1和很多无穷小ε组成。
那么其决定因素就为1,那么其行列式的值就为1了。
5、额外的问题与解答在之后的学习中发现一个问题就是我在猜想一步转移概率矩阵是否能收敛的问题上还是考虑的不够全面,漏掉了很多重要的问题,我也在这儿举例验证如下:Q=[0 1 0;0.5 0 0.5; 0 1 0]图(7):矩阵Q 的行列式值图(8):矩阵Q 的秩这个3阶的矩阵,也是书上的一个例题的矩阵,这个矩阵并不是上述我说的类单位矩阵或者是单位矩阵。
而是一个一般的矩阵,然而这个矩阵是没有办法收敛的其n 次的值是在两个值之间循环跳动的。
这个矩阵的det 值为0【见上图(7)】,但是并没有上述验证中的列相同达到收敛的规律。
但是其行列式的值也为0。
之后我算了一下他的秩,发现是2【见上图(8)】,也就是说秩的值小于阶的值,而我之前举的例子中,秩的值都是等于阶的值。
之后我又验证了一个矩阵【见下图(9)&图(10)】:W=[0.1 0.1 0.1 0.7;0 0.2 0.2 0.6;0 0 0.4 0.6;0.1 0.1 0.1 0.7]图(9):矩阵W 的行列式值图(10):矩阵W 的秩这是一个非满秩的矩阵所以他的行列式的值一定为0。
与我上述的结论冲突了,所以我上述的结论应建立在给出的一步转移概率矩阵为满秩的情况下才能成立。
若不为满秩的话,则可以算其各列的方差的平均值来进行比较,单位矩阵的列平均方差为(n-1)/n 而其他的一步转移概率矩阵则介于0-(n-1)/n 之间。
参考文献[1].曾勇等. <排队现象的建模、解析与模拟>.西安电子科技大学出版社.2011年9月 [2].对于一步转移矩阵收敛快慢的解答.豆丁网[3].吴广艳等.<MATLAB 简明实例教程>.东南大学出版社.2016年1月 [4]. Angle Roh 等.<转移概率矩阵>.MBA 智库网站.2009年 [5].居余马等.<线性代数第二版>.清华大学出版社.2002年9月作者简介:XXX :计算机学院计算机科学与技术专业,学号XXXXXXXXXXX ;XXX :计算机学院计算机科学与技术专业,学号 XXXXXXXXXXX第二题:分析多服务窗等待制M/M/N 排队系统,其中平均到达速率为l ,每个服务员的平均服务速率为μ由概率分布求系统中总顾客数的均值L . 虑到公式推导的复杂性,请用自己熟悉的语言“纸上写代码”,给出求解L 近似值的核心代码。
代码关键部分必须有注释.1. M /M /C 模型中,第一个M 表示顾客的到达为泊松流,第二个M 表示服务为独立同负指数分布,C 表示C 个服务员,系统容量为无穷,默认顾客源为无穷,排队规则采用FIFO(先到先服务)规则。
2. 令时刻t 系统内的总顾客数为N(t),取足够小的时间间隔,则单位时间N(t)只能加减一或者不变,则N(t)为生灭过程,取值范围为0到无穷 状态流图3.达到稳态时系统满足如下方程(列表也可)()()011111(1)()1()c c n n c n c n n c n P P P n P n P n c P cP C P n c μλμλμλμλμλ-+-+=++=+<=<+=+>=推导得:1()()!1()()!nc n n c n cn c n P n c C C λμλμ-⎧<⎪⎪=⎨⎪>=⎪⎩由归一性得:110011()()!!(1)C n c c c n P n C c λλρρμ--=⎡⎤=+⎢⎥-⎣⎦∑ .4. 编写代码求1110111()()(1)!(1)!(1)C n c c c n c L P n C c λλλμμρμ---=⎡⎤=+⎢⎥---⎣⎦∑的近似值,要体现精度控制。
代码如下(Java 语句): package paiduil;import java.math.*;import java.util.Scanner;public class test {public static void main(){System.out.println("请给定输入速率p:");Scanner in1=new Scanner(System.in);//输入数据和调用函数代码int p=in1.nextInt();System.out.println("请给定服务速率u:");Scanner in2=new Scanner(System.in);int u=in2.nextInt();System.out.println("请给定服务窗口数N:");Scanner in3=new Scanner(System.in);int N=in3.nextInt();System.out.println("solution(p,u,N)");}public static double x=0;double sum(double[] a)//数组求和{for(double e:a){x=x+e;}return x;}double factorial(int n)//阶乘{int i;int s=1;for(i=1;i<=n;i++)s*=i;return s*1.0;}double solution(double lambda,double u,int N){double p = lambda / u;double p_slash = lambda / (N * u);double[] factor_mat=new double[100];double[] P=new double[100];double length=0;int i;if (p_slash>=1){System.out.println("The result is infinite.\n");return -1;}for( i = 1 ;i<=(N+1);i++)//p[i]的系数factor_mat[i] = 1 / factorial(i-1) * (Math.pow(p,i-1));P[1] = 1 / (sum(factor_mat) + p_slash / (1-p_slash) *Math.pow(p , N)/factorial(N));//P[0]System.out.printf("P[0] = %.5f\n", P[1]);for( i = 2 ;i<= (N+1) ;i+=1)//p[i]P[i] = factor_mat[i] * P[1];System.out.printf("P[N]= %.5f\n", P[N+1]);//输出p[N]for( i = 1 ;i<=(N+1);i++)//0——N 队长length+=(i-1)* P[i];i = 1;while (true){double delta;double precision;//精度delta = (N + i) * Math.pow(p_slash, i) * P[N + 1];//P[N+1]即P[N]precision = P[N+1] * Math.pow(p_slash, i) * (N / (1-p_slash) + (i * (1-p_slash) + p_slash) / Math.pow((1-p_slash) , 2));if (precision < 1e-6)//达到精度{System.out.printf("已达到精度:%.8f\n", precision);System.out.println("此时队长:");break;}length = length + delta;System.out.printf("num %d: length = %.8f\n", i, length);//输出每一次循环判断精度的结果i=i+1;}return length;}}运行结果:。