发明名称 一种基于GSPN可靠性模型的可靠度分析方法
摘要 本发明公开了一种基于GSPN可靠性模型的可靠度分析方法,该方法有五大步骤:一、从GSPN模型的初态开始分析,将初始显态记录在可达矩阵ReachStaMat中;二、依次分析可达矩阵ReachStaMat中的每一状态,将每一状态可转移到的下一显态(集)记录到ReachStaMat中;三、由获得的矩阵GeneratMat求解连续时间马尔克夫链的转移速率矩阵;四、根据马尔克夫动态过程,由马尔克夫链的转移速率矩阵求解系统可靠度;五、分析结束。本发明是针对飞控系统GSPN可靠性模型的计算机自动化求解过程而提出的,这种方法节约了计算机的内存存储空间,提高了计算效率,同时也利于指出系统建模不当的地方。它在可靠性分析技术领域里具有实用价值和广阔的应用前景。
申请公布号 CN101846978A 申请公布日期 2010.09.29
申请号 CN201010184696.X 申请日期 2010.05.20
申请人 北京航空航天大学 发明人 孙晓哲;李卫琪;陈宗基;董卓宁
分类号 G05B13/04(2006.01)I 主分类号 G05B13/04(2006.01)I
代理机构 北京慧泉知识产权代理有限公司 11232 代理人 王顺荣;唐爱华
主权项 1.一种基于GSPN可靠性模型的可靠度分析方法,其实施对GSPN模型有以下要求:(1)GSPN模型的可达状态集有限;(2)GSPN模型中赋时变迁的触发概率服从指数分布函数<img file="FSA00000132875800011.GIF" wi="401" he="60" />分布函数中的实参数λ<sub>t</sub>即是飞控系统元件的故障率,因飞控系统可靠性要求高,其数量级一般为10<sup>-3</sup>以下;(3)在计算机上使用Matlab软件实现,并要求计算机的基本配置为CPU 2.80GHz、内存1G;其特征在于:该方法具体步骤如下:步骤一:从GSPN模型的初态开始分析,将初始显态记录在可达矩阵ReachStaMat中;显态指只使能赋时变迁、不使能瞬态变迁的状态;隐态为使能包含瞬态变迁的状态;步骤二:依次分析可达矩阵ReachStaMat中的每一状态,将每一状态可转移到的下一显态集记录到ReachStaMat中,并通过隐态消去的方法直接获得此状态转移到下一显态集的转移速率并存储到转移速率矩阵GeneratMat的相应位置;定义id1为ReachStaMat中状态索引号,令id1=1,步骤二1:分析ReachStaMa t的第id1个状态,记录id1显态下各使能变迁分别触发后转移到的状态及状态转移概率到TempStaMat和TempStaProb;依次分析TempStaMat中的每个状态,根据状态显隐态的判断,将显态记录到可达矩阵ReachStaMat,隐态记录在VashTransMat矩阵中进行再次触发分析;同时记录id1态向它可转移到的显态的转移速率到GeneratMat中;设GSPN的扩展可达集R中,T为显态集,V为隐态集,i,j表示扩展可达集中的任意两显态,即i,j∈T,从显态i向显态j的转移包括①从i到j的一步转移概率;②从显态i沿着一条全部由隐态构成的中间状态的路径,从隐态r转移到显态j;则显态i向显态j的转移概率u<sub>ij</sub>计算如下:<maths num="0001"><![CDATA[<math><mrow><msub><mi>u</mi><mi>ij</mi></msub><mo>=</mo><msub><mi>f</mi><mi>ij</mi></msub><mo>+</mo><munder><mi>&Sigma;</mi><mrow><mi>r</mi><mo>&Element;</mo><mi>V</mi></mrow></munder><msub><mi>P</mi><mi>r</mi></msub><mo>{</mo><mi>i</mi><mover><mo>&RightArrow;</mo><mi>r</mi></mover><mi>j</mi><mo>}</mo><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>1</mn><mo>)</mo></mrow></mrow></math>]]></maths>其中:f<sub>ij</sub>表示从显态i到j的一步转移概率,<img file="FSA00000132875800022.GIF" wi="191" he="118" />表示显态i沿着全部由隐态构成的中间状态,从隐态r首次转移到显态j的概率,其中路径可以包括任意步数;由转移概率u<sub>ij</sub>构造连续时间马尔克夫链的转移速率q<sub>ij</sub>,当i≠j,根据<maths num="0002"><![CDATA[<math><mrow><msub><mi>q</mi><mi>ij</mi></msub><mo>=</mo><munder><mi>lim</mi><mrow><mi>&Delta;t</mi><mo>&RightArrow;</mo><mn>0</mn></mrow></munder><mfrac><mrow><msub><mi>u</mi><mi>ij</mi></msub><mrow><mo>(</mo><mi>&Delta;t</mi><mo>)</mo></mrow></mrow><mi>&Delta;t</mi></mfrac><mo>,</mo></mrow></math>]]></maths>得:<maths num="0003"><![CDATA[<math><mrow><msub><mi>q</mi><mi>ij</mi></msub><mo>=</mo><munder><mi>lim</mi><mrow><mi>&Delta;t</mi><mo>&RightArrow;</mo><mn>0</mn></mrow></munder><mfrac><mrow><msub><mi>f</mi><mi>ij</mi></msub><mo>+</mo><munder><mi>&Sigma;</mi><mrow><mi>r</mi><mo>&Element;</mo><mi>V</mi></mrow></munder><msub><mi>P</mi><mi>r</mi></msub><mo>{</mo><mi>i</mi><mover><mo>&RightArrow;</mo><mi>r</mi></mover><mi>j</mi><mo>}</mo></mrow><mi>&Delta;t</mi></mfrac><mo>,</mo></mrow></math>]]></maths><maths num="0004"><![CDATA[<math><mrow><mo>=</mo><munder><mi>lim</mi><mrow><mi>&Delta;t</mi><mo>&RightArrow;</mo><mn>0</mn></mrow></munder><mfrac><msub><mi>f</mi><mi>ij</mi></msub><mi>&Delta;t</mi></mfrac><mo>+</mo><munder><mi>&Sigma;</mi><mrow><mi>r</mi><mo>&Element;</mo><mi>V</mi></mrow></munder><munder><mi>lim</mi><mrow><mi>&Delta;t</mi><mo>&RightArrow;</mo><mn>0</mn></mrow></munder><mfrac><mrow><msub><mi>P</mi><mi>r</mi></msub><mo>{</mo><mi>i</mi><mover><mrow><mo>&RightArrow;</mo><mo>j</mo></mrow><mtext>r</mtext></mover><mo>}</mo></mrow><mi>&Delta;t</mi></mfrac><mo>=</mo><msub><mi>q</mi><mrow><mi>ij</mi><mrow><mo>(</mo><mn>1</mn><mo>)</mo></mrow></mrow></msub><mo>+</mo><munder><mi>&Sigma;</mi><mrow><mi>r</mi><mo>&Element;</mo><mi>V</mi></mrow></munder><msub><mi>q</mi><mi>ij</mi></msub><mrow><mo>(</mo><mi>r</mi><mo>)</mo></mrow><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>2</mn><mo>)</mo></mrow></mrow></math>]]></maths>其中:q<sub>ij</sub>为由显态i向显态j的转移速率,q<sub>ij(1)</sub>为显态i到显态j的一步转移速率,q<sub>ij(r)</sub>为显态i经过隐态r转移到显态j的转移速率;连续时间马尔克夫链的状态转移速率矩阵Q是以q<sub>ij</sub>为元素的矩阵;分析ReachStaMat的第id1个状态可首次转移到的显态以及获得转移速率的具体步骤流程如下:定义id2为TempStaMat中状态索引号,令id2=1,步骤二1.1:分析TempStaMat中的第id2个状态,判断id2态下使能变迁是否全部为瞬态变迁,如果是,则此状态为隐态,转至步骤二1.3执行;否则为显态,转至步骤二1.2;步骤二1.2:状态id2为显态,记录此显态到ReachStaMat和id1到id2的一步转移速率到GeneratMat;将此显态添加在可达阵ReachStaMat前判断此显态是否已经记录在可达阵ReachStaMat中,若已记录,则不需要再次添加;显态id1转移到此显态id2的转移概率为TempStaProb<sub>id2</sub>,根据式(2),由<img file="FSA00000132875800031.GIF" wi="541" he="120" />获得显态id1到显态id2的一步转移速率q<sub>ij(1)</sub>;将q<sub>ij(1)</sub>添加到GeneratMat之前判断GeneratMat的相应位置是否已经存储了转移速率,若没有存储,直接添加;若已经存储,设q<sub>ij(原来)</sub>表示存储位置上原来存储的转移速率,执行q<sub>ij(原来)</sub>+q<sub>ij(1)</sub>,并将结果存储到相应位置;转至步骤二1.4;步骤二1.3:状态id2为隐态,将此隐态记录在隐态转移矩阵VashTransMat中的第一列,分析系统经过id2隐态可首次转移到的显态;定义m为VashTransMat中状态索引号,令m=1,步骤二1.3.1:分析VashTransMat的第m状态,若此状态为显态,则不进行分析,转至步骤二1.3.2;若此状态为隐态,根据隐态m下使能的瞬态变迁的优先级设置,触发优先级最高的瞬态变迁,并将触发后的系统状态进行是否存储判断后添加到VashTransMat中;同时记录m状态到此状态的转移概率到VashTransProb中的相应位置,矩阵VashTransProb是与VashTransMat中状态对应的状态转移概率矩阵;步骤二1.3.2:判断是否已经对VashTransMat中所有状态进行了分析,如果否,执行m=m+1,分析VashTransMat的下一个状态,转至步骤二1.3.1;如果是,顺序执行;步骤二1.3.3:分析转移概率矩阵VashTransProb,判断是否存在隐态吸收环;从隐态id2开始的可达状态集为VashTransMat,将可达集分为两个集合M<sub>T</sub>显态和M<sub>V</sub>隐态,对所有的状态进行重新排列,所有隐状态在前即隐态id2依旧放在第一列,显状态在后,转移概率矩阵VashTransProb根据可达集各状态位置之间的变换也相应的进行了重新排列,位置排列之后的转移概率矩阵为:<maths num="0005"><![CDATA[<math><mrow><msub><mi>P</mi><mi>vash</mi></msub><mo>=</mo><mfenced open='[' close=']'><mtable><mtr><mtd><msubsup><mi>P</mi><mi>v</mi><mi>VV</mi></msubsup></mtd><mtd><msubsup><mi>P</mi><mi>v</mi><mi>VT</mi></msubsup></mtd></mtr></mtable></mfenced></mrow></math>]]></maths>其中<img file="FSA00000132875800042.GIF" wi="75" he="58" />和<img file="FSA00000132875800043.GIF" wi="75" he="57" />分别表示矩阵VashTransMat包含的状态中隐态到隐态和隐态到显态的转移概率;设I表示与<img file="FSA00000132875800044.GIF" wi="76" he="58" />同维数的单位矩阵,判断矩阵<img file="FSA00000132875800045.GIF" wi="146" he="57" />是否奇异,若矩阵奇异,则从隐态id2开始触发之后系统存在隐态吸收环,建模错误,算法退出;转至步骤五;若矩阵非奇异,则系统不存在隐态吸收环;将可达显态M<sub>T</sub>添加到ReachStaMat矩阵之前,首先判断ReachStaMat中是否已经记录了显态M<sub>T</sub>,若没有记录,则直接将显态M<sub>T</sub>添加到ReachStaMat;以下计算显态id1经过隐态id2转移到显态集M<sub>T</sub>的转移概率U<sub>vash</sub>为:<maths num="0006"><![CDATA[<math><mrow><msub><mi>U</mi><mi>vash</mi></msub><mo>=</mo><msup><mrow><mo>(</mo><mi>I</mi><mo>-</mo><msubsup><mi>P</mi><mi>v</mi><mi>VV</mi></msubsup><mo>)</mo></mrow><mrow><mo>-</mo><mn>1</mn></mrow></msup><msubsup><mi>P</mi><mi>v</mi><mi>VT</mi></msubsup></mrow></math>]]></maths>转移概率矩阵U<sub>vash</sub>中与隐态id2对应的行向量即第一行U<sub>vash(1,:)</sub>为隐态id2转移到VashTransMat中显态集M<sub>T</sub>的转移概率;行向量TempStaProb<sub>id2</sub>*U<sub>vash(1,:)</sub>为显态id1经过隐态id2转移到显态M<sub>T</sub>的转移概率;根据式(2),由<img file="FSA00000132875800052.GIF" wi="726" he="116" />获得显态id1经过隐态id2转移到显态集M<sub>T</sub>的转移速率q<sub>ij(r)</sub>;将q<sub>ij(r)</sub>添加到GeneratMat之前判断GeneratMat的相应位置是否已经存储了转移速率,若没有存储,直接添加;若已经存储,则执行q<sub>ij(原来)</sub>+q<sub>ij(r)</sub>,并存储到相应位置;步骤二1.3.4:清空临时矩阵VashTransMat和VashTransProb,节约存储空间;步骤二1.4:判断是否对临时矩阵TempStaMat中的状态全部分析完成,若是,则清空临时矩阵TempStaMat和TempStaProb,顺序执行;若否,则执行id2=id2+1,转至步骤二1.1执行;步骤二2:判断是否分析完ReachStaMat中的所有状态,若是,则顺序执行;若否,执行id1=id1+1,转至步骤二1执行;步骤三:由至此获得的矩阵GeneratMat求解连续时间马尔克夫链的转移速率矩阵;连续时间马尔克夫链的转移速率矩阵中,每行元素之和均为零,即<img file="FSA00000132875800053.GIF" wi="208" he="127" />l为马尔克夫过程的状态数目;至此,经过以上算法后存储的矩阵GeneratMat中只包含了不同状态之间的转移速率,即非对角线元素;根据<img file="FSA00000132875800061.GIF" wi="325" he="128" />由矩阵GeneratMat求解状态转移到自身的转移速率;步骤四:根据马尔克夫动态过程,由马尔克夫链的转移速率矩阵求解系统可靠度;设可达阵ReachStaMat的列数为l,概率向量P(t)=(p<sub>1</sub>(t),p<sub>2</sub>(t),...,p<sub>l</sub>(t)),其中p<sub>i</sub>(t)为系统处于状态T<sub>i</sub>的瞬态概率,则有下面的微分方程成立:<maths num="0007"><![CDATA[<math><mrow><mfenced open='{' close=''><mtable><mtr><mtd><mrow><mo>(</mo><msub><msup><mi>p</mi><mo>&prime;</mo></msup><mn>1</mn></msub><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow><msub><msup><mi>p</mi><mo>&prime;</mo></msup><mn>2</mn></msub><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow><mo>,</mo><mo>.</mo><mo>.</mo><mo>.</mo><mo>,</mo><msub><msup><mi>p</mi><mo>&prime;</mo></msup><mi>l</mi></msub><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow><mo>)</mo></mrow><mo>=</mo><mrow><mo>(</mo><msub><mi>p</mi><mn>1</mn></msub><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow><mo>,</mo><msub><mi>p</mi><mn>2</mn></msub><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow><mo>,</mo><mo>.</mo><mo>.</mo><mo>.</mo><mo>,</mo><msub><mi>p</mi><mi>l</mi></msub><mrow><mo>(</mo><mi>t</mi><mo>)</mo></mrow><mo>)</mo></mrow><mo>*</mo><mi>GeneratMat</mi></mtd></mtr><mtr><mtd><mi>P</mi><mrow><mo>(</mo><mn>0</mn><mo>)</mo></mrow><mo>=</mo><mo>[</mo><msub><mi>p</mi><mn>1</mn></msub><mrow><mo>(</mo><mn>0</mn><mo>)</mo></mrow><mo>,</mo><msub><mi>p</mi><mn>2</mn></msub><mrow><mo>(</mo><mn>0</mn><mo>)</mo></mrow><mo>,</mo><mo>.</mo><mo>.</mo><mo>.</mo><mo>,</mo><msub><mi>p</mi><mi>l</mi></msub><mrow><mo>(</mo><mn>0</mn><mo>)</mo></mrow><mo>]</mo></mtd></mtr></mtable></mfenced><mo>-</mo><mo>-</mo><mo>-</mo><mrow><mo>(</mo><mn>3</mn><mo>)</mo></mrow></mrow></math>]]></maths>系统在初始时刻的状态标识为T<sub>i</sub>的概率为p<sub>i</sub>(0),由系统的结构以及初始时刻的条件确定P(0);求解上面的微分方程,获得每个显态的瞬态概率;由系统的失效状态集合求解系统失效概率,则得系统可靠度;步骤五:分析结束。
地址 100191 北京市海淀区学院路37号