CN103440393B - 一种面向连续时间马尔科夫链的状态空间约简方法 - Google Patents

Abstract

本发明公开了一种面向连续时间马尔科夫链的状态空间约简方法，包括以下步骤：在局部状态空间上配置连续随机逻辑中X，U，R，P，W算子的限界语义；将连续时间马尔科夫链转换为离散时间马尔科夫链，并利用均匀化技术计算局部空间上的瞬态与稳态概率；对于算子X，U，R将时域分成零点到某个时间点，两个时间点之间，某个时间点到无穷大三个时间段，并分别利用瞬态概率计算算子P对应的概率度量；对于稳定算子W，构造一组线性不等式来约束离散时间马尔科夫链中瞬态概率与稳态概率之间的关系，求出该不等式的解作为稳态概率的度量值。本发明只需遍历分析属性所需的局部空间，可有效约简状态空间，可应用于大规模随机系统的性能与可靠性分析。

Description

一种面向连续时间马尔科夫链的状态空间约简方法

技术领域

本发明属于随机系统性能与可靠性分析技术领域，涉及面向连续时间马尔科夫链的状态空间约简技术。

背景技术

模型检测是一种自动化程度非常高的有限状态系统验证技术，目前已经在计算机硬件、通信与安全协议、软件可靠性的验证方面获得了较大的成功。传统模型检测技术关注的是系统行为的绝对正确性，如系统不能进入死锁状态。然而分布式算法，多媒体协议，容错系统等往往关心某种量化属性，如消息传送失败的概率不高于1%，在时间t内至多m个消息丢失的概率不高于0.8%，请求发送后在5到7个时间单元内得到响应的概率不低于70%等等。随机模型检测致力于解决这类属性的自动化验证问题。

在随机模型检测中一般使用概率计算树逻辑PCTL和连续随机逻辑CSL刻画属性，使用马尔科夫过程建立系统模型，主要包括离散时间马尔可夫链，马尔科夫决策过程，连续时间马尔科夫链等。每种模型都具有一定的特性，不同的特性决定了模型表达和分析的重点不一样。与其它模型相比，连续时间马尔科夫链的主要特性在于能刻画连续时间和指数分布。这两种特性使得连续时间马尔科夫链的模型检测近年来成为一种成功的定量分析技术。

连续时间马尔科夫链的模型检测技术主要关注于随机系统的性能、可靠性等性质的定量分析。例如Ender Yüksel通过利用连续时间马尔科夫链为我国智能电网中的传感网络建立随机模型，并计算出了长期运行中传感节点失效的概率，以及更换传感器中电池的最优时间段。Shinji Kikuchi对云计算系统中并发实时迁移操作的性能进行了分析，计算出了某个时间段内发送服务器上多于4个迁移操作的概率，以及某个时间段内接受服务器上超过3个迁移操作正在处理的概率。M.Kwiatkowska通过分别计算传感器、执行器、输入输出处理器、中心计算处理器引起系统关闭的概率，分析了嵌入式控制系统的可靠性。

连续时间马尔科夫链的模型检测技术在生物学领域也有着重要的应用。MartaZ.Kwiatkowska分析了成纤维细胞生长因子信号通路的健壮性，并给出了系统各种动态行为的量化度量，从而加深了对信号通路的理解。J.Heath对分裂素激活的蛋白激酶级联反应系统中各个成分之间的交互进行了定量刻画。这些成功的应用实例说明模型检测连续时间马尔科夫链是对马尔科夫过程传统分析技术的有力扩展与补充。

目前连续时间马尔科夫链的模型检测方法是一种全局检测方法，即通过遍历整个系统的全局空间完成属性的分析，因此与传统模型检测一样，状态空间爆炸依旧是模型检测连续时间马尔科夫链实用化的主要瓶颈（这里状态空间爆炸是指对于并发系统，其状态的数目往往随着并发分量的增加呈指数增长），约简状态空间对提高模型检测连续时间马尔科夫链技术的实用性至关重要。

限界模型检测是一种有效的空间约简方法，其基本思想是在有限的局部空间中逐步搜索属性成立的证据或者失效的反例，从而达到约简状态空间的目的。连续时间以及指数分布是连续时间马尔科夫链上的两个主要特性，这两种特性为应用限界检测来约简连续时间马尔科夫链上的状态空间带来了新的重要的问题，如在有穷路径约束下瞬态概率与稳态概率的计算，不同的时间约束类型对限界满足性的影响不同，以设置路径长度的上限作为判断算法终止的标准已经失效，必须设计新的标准等等。这些新问题说明将限界检测方法应用于约简连续时间马尔科夫链是不平凡的，有必要对该方法进行全新的研究。

发明内容

本发明的目的在于提供一种面向连续时间马尔科夫链的状态空间约简方法，以提高模型检测连续时间马尔科夫链在分析随机系统性能与可靠性上的实用性，提高可处理系统的规模。

为了解决以上技术问题，本发明的采用的技术方案如下。

一种面向连续时间马尔科夫链的状态空间约简方法，其特征在于包括以下步骤：

步骤一，在局部状态空间上配置连续随机逻辑中X，U，R，P，W算子的限界语义

令为连续时间马尔科夫链。对任意的状态s∈S和界k，φ∈CSL，满足性关系s|＝_kφ递归配置如下：

·对任意的s∈S，s|＝_ktrue；

·s|＝_ka当且仅当a∈L(s)；

·当且仅当

·当且仅当s|＝_kφ，s|＝_kψ；

·当且仅当s|＝_kφ或者s|＝_kψ；

·s|＝_kP_～p[φ]当且仅当Prob^C(s,φ,k)～p；

·s|＝_kW_～p[φ]当且仅当 ${\mathrm{\Σ}}_{s^{\′}|=\mathrm{\φ}}{\mathrm{\π}}_s^C(s^{\′},k)~p,$

其中Prob^C(s,φ,k)＝Pr_s{ω∈Path^C(s)|ω|＝_kφ}，且对任意的路径ω∈Path^C(s)：

ω|＝_kXφ当且仅当k≥1，且ω(1)|＝_kφ；

ω|＝_kφU^Iψ当且仅当

ω|＝_kφR^Iψ当且仅当或者

步骤二，限界下计算局部空间上状态之间的k界瞬态概率

连续时间马尔科夫链上的嵌入离散时间马尔科夫链配置为其中对S中任意的状态s，s'：

连续时间马尔科夫链上的矩阵配置如下；

对任意的连续时间马尔科夫链和矩阵Q，均匀化的离散时间马尔科夫链配置为其中P^unif(C)＝I+Q/q，q是比max{E(s)|s∈S}大的任意整数。

k界瞬态概率向量计算方法如下： ${\mathrm{\Π}}_{t,k}^C=\underset{i=0}{\overset{k}{\mathrm{\Σ}}}{\mathrm{\γ}}_{i,q,t}\·{\left(P^{\mathrm{unif}\left(C\right)}\right)}^i,$ 其中 ${\mathrm{\γ}}_{i,q,t}=e^{-q\·t}\·\frac{{(q\·t)}^i}{i!}.$ 该计算方法直观上的解释是：均匀化的离散时间马尔科夫链中每一步都对应着参数为q的指数分布延迟，矩阵的幂运算(P^unif(C))ⁱ给出了在离散时间马尔科夫链中在i步内状态之间的转换概率，γ_i,q,t是参数为q·t的泊松过程。

步骤三，限界下计算局部空间上X，U，R，P，W算子的概率度量

输入：连续时间马尔科夫链CSL公式β，以及界k。

输出：在界k的约束下所有满足β的状态的集合Sat(β)＝{s∈S|s|＝_kβ}。

计算流程：首先递归地计算满足每个子公式的状态的集合，最后得出集合Sat(β)。Sat(β)的计算过程如下：

Sat(true)＝S;

Sat(a)＝{s|a∈L(s)}；

Sat(P_～p[φ])＝{s∈S|Prob^C(s,φ,k)～p}；

$\mathrm{Sat}\left(S_{~p}\right[\mathrm{\φ}\left]\right)=\{s\∈S|{\mathrm{\Σ}}_{s^{\′|=\mathrm{\φ}}}{\mathrm{\π}}_s^C(s^{\′},k)~p\}.$

上述计算过程除了P_～p[·]和S_～p[·]算子，其它算子的处理是平凡的，因此现在给出这两个算子当中概率度量的计算。

情形1公式β形式为P_～p[Xφ]，且集合Sat(φ)已知

X算子与连续时间马尔科夫链中的时间延迟无关，仅仅依赖于从当前状态转移到下一个状态的概率，因此概率度量Prob^C(s,Xφ,k)计算如下：

情形2公式β形式为P_～p[φU^Iψ]，且集合Sat(φ)，Sat(ψ)已知

对于该算子我们需要对每个状态s计算概率Prob^C(s,φU^Iψ,k)，其中I是任意的非负实数区间。注意到Prob^C(s,φU^Iψ,k)＝Prob^C(s,φU^cl(I)ψ,k)，其中cl(I)是区间I的闭包，且Prob^C(s,φU^[0,∞)ψ,k)＝Prob^emb(C)(s,φU^≤∞ψ,k)。具体分下列三种情况讨论Prob^C(s,φU^Iψ,k)的计算：

·I＝[0,t]，这里

·I＝[t,t']，这里且t≤t'；

·I＝[t,∞)，这里

情形2.1I＝[0,t]

对任意的连续时间马尔科夫链和CSL公式φ，令为连续时间马尔科夫链，其中对任意不满足φ的状态s，R[φ](s,s')＝R(s,s')，否则R[φ](s,s')＝0。

对于连续时间马尔科夫链和CSL公式φ，ψ，以及正实数Prob^C(s,φU^[0,t]ψ,k)的计算方法为

直观上上述计算方法可以解释为在连续时间马尔科夫链中，一旦到达一个满足ψ的状态，路径就不能退出这个状态，且一旦进入满足的状态，就永远不能到达满足ψ的状态，所以路径公式φU^[0,t]ψ在k界的约束下被满足的概率等价于在连续时间马尔科夫链中在k界的约束下在时刻t处于满足φ的状态概率。

情形2.2I＝[t,t']

依据ω|=_kφU^[t,t']ψ的语义解释，在这种情况下可以将路径ω分解成两个部分：(a)在时间t内处于满足φ的状态；(b)在时间t'-t内在到达满足ψ的状态，之前一直处于满足φ的状态。对于前者，我们利用I＝[0,t]当中的思想通过计算连续时间马尔科夫链中的转移概率得到相应的概率度量。Prob^C(s,φU^[t,t']ψ,k)的计算如下：

情形2.3I＝[t,∞)

依据ω|=_kφU^[t,∞)ψ的语义解释，在这种情况下可以将路径ω分解成两个部分：(a)在时间t内一直处于满足φ的状态；(b)在之后的时间内到达满足ψ的状态，且到达之前一直处于满足φ的状态。第二个部分对时间没有约束，因此嵌入离散时间马尔科夫链在这种情况下可被使用。Prob^C(s,φU^[t,∞]ψ,k)的计算过程如下：

情形3公式β形式为P_～p[φR^Iψ]，且集合Sat(φ)，Sat(ψ)已知

情形3.1I＝[0,t]

依据ω|=_kφR^[0,t]ψ的语义解释，ω只需满足两种路径约束的一种即可，因此概率度量Prob^C(s,φR^[0,t]ψ,k)需分成两部分计算。我们注意到存在同时满足两种约束的路径，因此为了避免重复计算，下面定义一种结构。

对任意的连续时间马尔科夫链和CSL公式φ，令为连续时间马尔科夫链，其中对任意满足φ的状态s，R<φ>(s,s')＝R(s,s')，否则R[φ](s,s')＝0。

对于连续时间马尔科夫链和CSL公式φ，ψ，以及正实数

直观上Prob^C(s,φR^[0,t]ψ,k)的计算解释为：1）在连续时间马尔科夫链中，一旦到达一个满足φ的状态，路径就不能退出这个状态，且一旦进入满足的状态，就永远不能到达满足φ的状态；2）在连续时间马尔科夫链中，路径一直处于满足的状态。

情形3.2I＝[t,t']

依据ω|=_kφR^[t,t']ψ的语义解释，将路径ω分解成两个部分：1）在时间t内处于满足ψ的状态；2）在时间t'-t内到达满足φ的状态，之前（包括现在）一直处于满足ψ的状态，或者在[t,t']区间内一直处于满足ψ的状态。依据这种分解，概率度量Prob^C(s,φR^[t,t']ψ,k)计算如下：

情形3.3I＝[t,∞)

这种情况几乎等同于I＝[t,t']的情况。依据ω|=_kφR^[t,∞)ψ的语义解释，将路径ω分解成两个部分：1）在时间t内处于满足ψ的状态；2）在此后的时间内到达满足φ的状态，到达之前（包括到达时刻）一直处于满足ψ的状态，或者在[t,∞)区间的一直处于满足ψ的状态。

概率度量Prob^C(s,φR^[t,∞)ψ,k)计算如下：

情形4公式β形式为W_～p[φ]，且集合Sat(φ)已知

令为连续时间马尔科夫链，为嵌入离散时间马尔科夫链。对于状态s∈S,1)如果s₁＝s，则称s₁是从s出发0步可达的；2)如果s_l-1是从s出发l-1步可达的，且P^emb(C)(s_l-1,s_l)＞0，则称s_l是从s出发l步可达的。

引入记号S^s,k表示从状态s出发k步内可达的状态组成的集合。对于引入记号表示s的前驱状态的集合。记号x(s,k)表示

对任意的 $s\&Element;S^{\stackrel{\&OverBar;}{s},k},$ 令 $x(s,k)\&le;\underset{s^{\&prime;}\&Element;\mathrm{pre}\left(s\right)}{\mathrm{\&Sigma;}}x\left(s^{\&prime;}\right)\&CenterDot;P^{\mathrm{emb}\left(C\right)}(s^{\&prime;},s)---\left(1\right).$

对于上述不等式（1），求出该不等式的最小解，从而得到稳态概率的下近似值。

各参数的解释如下：

本发明具有有益效果。模型检测连续时间马尔科夫链技术可以用来分析随机系统的性能与可靠性，状态空间爆炸是应用该技术处理大规模随机系统的主要瓶颈。本发明提出一种在有限的局部空间中逐步搜索属性成立的证据的限界检测技术，其特点在于只遍历分析属性所需的局部空间，因此能够有效克服状态空间爆炸，显著提高可分析随机系统的规模。

具体实施方式

下面结合实施例对本发明的技术方案做进一步详细说明。

锅炉控制系统组成部分包括：主处理器M，输入处理器I，输出处理器O，三个温度传感器S1，S2，S3，两个阀门A1，A2。处理器之间通过总线连接，传感器采集数据后将数据发送给输入处理器I融合，融合完毕后I将数据发送给M，M计算后将命令发送给处理器O，O再依据命令控制阀门A1，A2。

任何温度传感器都可能失效，对于输入处理器而言最多允许一个温度传感器失效，当多于一个温度传感器失效时，输入处理器通知主处理器关闭系统。执行器也可能失效，对于输出处理器而言最多允许一个执行器失效，当多于一个执行器失效时，输出处理器通知主处理器关闭系统。输入输出处理器本身也会失效，这种失效可能是暂时的，也可能是永久的。无论哪种失效，都会导致主处理器M无法从I读取数据，无法发送指令给O，从而使得M跳出当前的执行周期。如果M连续跳出执行周期的次数超过MAX_COUNT，则M会关闭系统。

PRISM是由牛津大学的Marta Kwiatkowska教授主导开发，一款面向学术界可免费使用的概率模型检测工具，主要用来对随机系统的行为进行建模与分析。我们采用PRISM中的建模语言为锅炉控制系统建立其连续时间马尔科夫链模型，并验证相关属性。

考察温度传感器引起系统关闭的可能性，具体验证属性：温度传感器引起系统关闭的概率不低于0.5。该属性利用概率模型检测工具PRISM中的属性表示方式可以描述为：P＞＝0.5[!downU^[0,^∞]fail_sensors]，这里down表示系统关闭，fail_sensors表示温度传感器失效。

令k为局部空间的深度。

第一步，依据上述描述的状态之间的转换关系，构造可达深度为k的连续时间马尔科夫链C；

第二步，构造连续时间马尔科夫链C[P＞＝0.5[!downU^[0,∞]fail_sensors]]；

第三步，构造嵌入离散时间马尔科夫链emb(C)；

第四步，对于初始状态s₀计算概率度量Prob^C(s₀,!downU^[0,∞]fail_sensors,k)。

在不同MAX_COUNT取值情况下本发明与全局检测算法所需空间的对比情况如表1所示。

表1.在不同MAX_COUNT取值情况下本发明与全局检测算法所需空间的对比

从表中不难发现本发明只需遍历一半的状态空间即可完成可靠性分析，因此本发明在确保完成属性分析的前提下有效约简了状态空间，提高了可分析系统的规模，可应用于大规模随机系统的分析。

Claims (1)

1.一种面向连续时间马尔科夫链的状态空间约简方法，其特征在于包括以下步骤：

步骤一，在将路径长度约束在有限长度的情况下配置连续随机逻辑中X，U，R，P，W五个算子的限界语义；X为表示“下一个时刻”的时态算子，U为表示“一直到”的时态算子，R为表示“释放”的时态算子，P为概率度量算子，W为稳态算子；具体配置按照以下步骤执行：

步骤1,配置P算子的限界语义：s|＝_k P_～p[φ]当且仅当Prob^C(s,φ,k)～p；

步骤2,配置W算子的限界语义：s|＝_k W_～p[φ]当且仅当

步骤3,计算概率度量Prob^C(s,φ,k)：Prob^C(s,φ,k)＝Pr_s{ω∈Path^C(s)|ω|＝_kφ}，且对任意的路径ω∈Path^C(s)：

ω|＝_k Xφ当且仅当k≥1，且ω(1)|＝_kφ；

ω|＝_kφU^Iψ当且仅当

ω|＝_kφR^Iψ当且仅当或者

步骤二，在将路径长度约束在有限长度的情况下计算局部空间上状态之间的k界瞬态概率，这里k是自然数；具体计算过程按照以下步骤执行：

步骤1,将连续时间马尔科夫链上的嵌入离散时间马尔科夫链配置为其中对S中任意的状态s，s'：

步骤2,将连续时间马尔科夫链上的矩阵配置为：

步骤3,对任意的连续时间马尔科夫链和矩阵Q，将均匀化的离散时间马尔科夫链配置为P^unif(C)＝I+Q/q，q是比max{E(s)|s∈S}大的任意整数；

步骤4,k界瞬态概率向量计算方法为

步骤三，在将路径长度约束在有限长度的情况下利用k界瞬态概率计算X，U，R，P，W算子对应的概率度量，计算过程具体为：

输入：连续时间马尔科夫链CSL公式β，以及界k；

输出：在界k的约束下所有满足β的状态的集合Sat(β)＝{s∈S|s|＝_kβ}；

计算流程：首先计算满足每个子公式的状态的集合，最后得出集合Sat(β)；Sat(β)的计算过程如下：

Sat(true)＝S；

Sat(a)＝{s|a∈L(s)}；

Sat(φ∧ψ)＝Sat(φ)∩Sat(ψ)；

Sat(φ∨ψ)＝Sat(φ)∪Sat(ψ)；

Sat(P_～p[φ])＝{s∈S|Prob^C(s,φ,k)～p}；

$Sat(W_{~p}\&lsqb;\mathrm{\&phi;}\&rsqb;)=\{s\&Element;S|{\mathrm{\&Sigma;}}_{s^{\&prime;}|=\mathrm{\&phi;}}{\mathrm{\&pi;}}_s^C(s^{\&prime;},k)~p\};$

各参数的解释如下：

φ、ψ是连续随机逻辑公式；

s表示连续时间马尔科夫链中的状态；

|＝_k是在界k的约束下状态满足连续随机逻辑公式；

～p表示p是[0,1]上的实数，～∈{＜,≤,＞,≥}；

Prob^C(s,φ,k)表示连续时间马尔科夫链C中在路径长度为k的约束下状态s满足连续随机逻辑公式φ的概率；

表示连续时间马尔科夫链C中在时刻t系统处于状态s的k界瞬态概率；

Pr_s{ω∈Path^C(s)|ω|＝_kφ}表示从状态s出发所有长度为k且满足φ的路径集合的概率；

ω表示一个非空序列s₀t₀s₁t₁s₂...，其中对任意的i≥0，R(s_i,s_i+1)＞0，t_i为非负实数；

ω(i)表示序列ω上的第i个状态；

Path^C(s)表示连续时间马尔科夫链C上从状态s出发的路径的集合，包括有穷路径和无穷路径；

I表示实数上的区间；

t表示时间；

ω@t表示对于无穷路径ω在时刻t系统所处的状态，即ω(j)，j是使得成立的最小的整数；

position(ω,t)表示对于路径ω＝s₀t₀s₁t₁s₂...，引入记号position(ω,t)表示使得成立的最小的整数j；表示“存在一个”；表示“任意的”；

表示连续时间马尔科夫链，其中S是有限状态集，是初始状态，是转移率矩阵(表示非负实数)，L:S→2^Ap是状态标记函数，为每个状态指定该状态下为真的原子命题；

表示嵌入离散时间马尔科夫链；

E(s)＝∑_s'∈SR(s,s')是从每个状态出发的转移率之和；

表示上的矩阵；

表示均匀化的离散时间马尔科夫链；

q表示比max{E(s)|s∈S}大的任意整数；

表示k界瞬态概率组成的向量；

γ_i,q,t表示参数为q·t的泊松过程；

Sat(β)＝{s∈S|s|＝_kβ}表示在界k的约束下所有满足β的状态的集合；

true表示永真的原子命题；

表示命题逻辑中的否定算子；

φ∧ψ表示公式φ，ψ之间的合取；

φ∨ψ表示公式φ，ψ之间的析取；

Sat(φ)∩Sat(ψ)表示集合Sat(φ)与Sat(ψ)之间的交集；

Sat(φ)∪Sat(ψ)表示集合Sat(φ)与Sat(ψ)之间的并集；

表示在连续时间马尔科夫链C中，从状态s出发，在长期的运行中处于状态s'的k界稳态概率；

t'表示时间；

s'，s"是连续时间马尔科夫链中的状态；

P^emb(C)表示嵌入离散时间马尔科夫链emb(C)上的转移率矩阵；

Ap表示原子命题的集合；

R(s,s')表示连续时间马尔科夫链中从状态s到状态s'的转移率；

R(s,s")表示连续时间马尔科夫链中从状态s到状态s"的转移率；

S\Sat(φ)表示从集合S中剔除Sat(φ)中的所有元素；

L(s)表示连续时间马尔科夫链中的函数L作用于状态s的值。