故障诊断技术对推进系统是至关重要的，但故障样本稀少且信息不全、诊断技术不够成熟，知识库维护困难等问题，使得诊断效率和准确性需要不断提高进步^[1-2]。国外诊断系统一般采用基于定性模型和多信号结合的诊断方法，根据建立的数学系统模型和输入信息，得到预测输出与实际输出对比，对系统进行故障诊断。这种方法不需要历史经验同时又能预测未知故障，但对模型精度要求很高，模型复杂^[3]。文献[4]提出建立故障传播的有向图，在有向图中加入故障触发概率信息和时间约束来对故障集合进一步刷选判断。文献[5-6]提出基于符号有向图（signed directed graph, SDG）模型和模糊推理相结合的故障诊断模型，利用故障概率信及其对应的触发权重对候选故障集进行故障可能性排序，提高了诊断精度。文献[7]利用贝叶斯网络推理算法和区间灰数运算规则给出不同故障下系统关键环节的触发概率和部件故障权重等评判标准，为不确定环境下的系统可靠性分析提供了理论依据。

基于定性模型的故障诊断方法存在不足，需要与其他方法结合来提高诊断方法的可行性，本研究提出部分可观时间Petri网的贝叶斯故障诊断方法^[8-22]，利用可观变迁时间序列信息和不可观节点触发的条件概率，计算故障节点触发概率，判断系统故障状态。

1 部分可观时间Petri网的贝叶斯故障检测 1.1 部分可观时间Petri网

定义1　Petri网（Petri nets，PN）为一个4元数组 ${{N}} = ({{S}},{{T}},{\rm{Pre}},{\rm{Post}})$ ， ${{S}}$ 是 $m$ 个库所（Place）的集合， ${{T}}$ 是 $n$ 个变迁(Transition)的集合， ${\rm{Pre}}:{{S}} \times {{T}} \to {{N}}$ 和 ${\rm{Post}}:{{S}} \times {{T}} \to {{N}}$ 分别是输入弧和输出弧的关联函数。部分可观时间Petri网（partially obesrved time petri nets, POTPN）定义为 ${{N}}_p = ({{N}},{{Q}},{{L}})$ ，其中 ${{N}} = ({{S}},{{T}},{\rm{Pre}},{\rm{Post}})$ 是Petri网的4元组， ${{Q}}:{{T}} \to $ ${{Q}} \times$ $ ({{Q}} \cup \{ \infty \} )$ 是变迁初始时间约束相关联的闭区间集合 ${{Q}}(t_i) = (l_i,u_i),l_i \geqslant 0,u_i \geqslant l_i$ ， $l_i \ne \infty $ ，当且仅当时间区间 $[l_i,u_i]$ 内,变迁 $t_i$ 才能触发。标签函数 $L:T \to L \cup$ $ \{ \varepsilon \} $ 为每个变迁分配一个标签 $L$ 或空符号 $\varepsilon $ 。不可观标签和静止标签集合为 $T_u$ , $T_u = \{ t \in T|{\cal{L}}(t) = \varepsilon \} $ ，可观标签集合为 $T_o = \{ t \in T|{\cal{L}}(t) = L\} $ 。

$\sigma = (t_0,\tau_ 0)(t_1,\tau_ 1)(t_2,\tau_ 2) \cdot \cdot \cdot (t_h,\tau_ h) \in (T \times {\bf{R}}_0^ + )$ 为变迁时间序列(transition time sequence, TTS),下标为 $i$ 的变迁 $t_i \in T$ 对应下标为 $i$ 的变迁触发时间点 $\tau_ i$ ，记为 $M_0[t_0(\tau_ 0)\rangle M_1[t_1(\tau_ 1)\rangle M_2 ,\cdot \cdot \cdot ,[t_h(\tau_ h)\rangle M_k,t_0 \leqslant t_1 \leqslant , \cdot \cdot \cdot , \leqslant t_h$ ，简记为 $M_0[\sigma \rangle M_k$ 或 $M_0[\sigma \rangle$ ，其中M₀为初始标识。

时间标签序列（time label sequence, TLS） $\delta_ o = (\gamma_{ e0},\tau_{ e0})(\gamma_{ e1},\tau_{ e1}), \cdot \cdot \cdot, (\gamma_{ eh},\tau_{ eh}) \in (L \times {\bf{R}}_0^ + )$ 是可观测变迁和对应触发时间的序列，任意 $\delta_ 0$ 至少存在一个变迁时间序列（TTS） $\sigma $ 与其对应。

1.2 部分可观时间贝叶斯Petri网

区别于一般Petri网，贝叶斯Petri网是针对故障变迁的有向无环网络，相关联节点之间都有触发的条件概率，同时贝叶斯Petri网可以省去部分与故障无关变迁和库所。贝叶斯网络 $B(G,P)$ 的组成包括n个节点的有向无环图 $G$ 和各节点对应的条件概率表 $P$ 。论域 ${{U}} = $ $\{ x_1,x_2, \cdots, x_n\} $ 为 $G$ 的节点变量， $p(x_i|{\rm{parents}}(x_i))$ 表示节点 $x_i$ 与父节点 ${\rm{parents}}(x_i)$ 中变量集合的条件概率。贝叶斯网络中全部变量集合的联合概率分布为

$p({{U}}) = p(x_1,x_2, \cdots ,x_n) = \prod\limits_{i = 1}^n {p(x_i|{\rm{parents}}(x_i))} $

(1)

对于随机变量 $x_i$ 的m 个基本事件 $\{ x_{i1},x_{i2}, \cdots, $ $x_{im}\} $ ，若获得的观察结果为 ${{V}} = \{ x_1,$ $x_2 ,\cdots, x_i - 1,x_i + 1, \cdots,$ $ x_n\} $ ，其条件概率为

$p(x_i = x_{ij}|V) = \frac{{p(x_i = x_{ij},{{V}})}}{{\displaystyle\sum\limits_{k = 1}^m {p(x_i = x_{ik},{{V}})} }}$

(2)

运用贝叶斯网络的逆向推理，由式(1)、(2) 可计算 $x_i = x_{ij}$ 的后验概率：

$\begin{array}{l} p(x_i = x_{ij}|x_1, \cdots ,x_i - 1,x_i + 1,x_n) = \\ \frac{\displaystyle{p(x_i = x_{ij})\prod\limits_{a \ne i} {p(x_a|{\rm{parents}}(x_a))} }}{{\displaystyle\sum\limits_{k = 1}^m {[p(x_i = x_{ik})\prod\limits_{a \ne i} {p(x_a|{\rm{parents}}(x_a))} ]} }} \end{array}$

(3)

定义2　贝叶斯Petri网(bayesian Petri net, BPN)定义为 ${{N}}_B = ({{S}},{{T}},{{P}},{{F}},{{B}},{{M}}_0)$ ，其中 ${{S}} = \{{{ S}}_N,{{S}}_B\} $ 是库所（Place）的有限集， ${{S}}_N = \{ s_{N_0},s_{N_1}, \cdots ,s_{N_m}\} $ 为常规库所有限集， ${{S}}_B = \{ s_{B_0},s_{B_1}, \cdots, s_{B_m}\} $ 为贝叶斯库所有限集； ${{T}} = \{ {{T}}_N,{{T}}_B\} $ 是变迁(Transition)的有限集， ${{T}}_N = $ $\{ t_{N_0},t_{N_1}, \cdots ,t_{N_m},\} $ 为常规变迁有限集， ${{T}}_B = \{ t_{B_0},t_{B_1}, \cdots,$ $ t_{B_m}\} $ 为贝叶斯变迁有限集； ${{F}} \subset ({{P}} \times {{T}}) \cup ({{T}} \times {{P}})$ 是输入弧和输出弧的有限集； ${{P}}:{{T}}_B \to [0,1],{{P}}(t_{B_m})$ 为贝叶斯变迁 $t_{B_m}$ 的条件概率表， ${{P}}(s_{B_m})$ 为贝叶斯库所 $s_{B_m}$ 的条件概率表； ${{B}}:{{T}}_B \to [0,1]$ ， ${{B}}(t_{B_m})$ 为 ${}^ \bullet t_{B_m} $ 在 ${}^ \bullet t_{B_m}$ 条件下的概率分布函数， $t_{B_m}^ \bullet $ 为 $t_{B_m}$ 的后置库所集。

$\begin{array}{l} {{B}}(t_{B_i}) = [b(s_1),b(s_2), \cdots, b(s_j)] \\ b(s_i) = p(si_|{}^ \bullet t_{B_i}),i = 1,2, \cdots ,j \end{array} $

(4)

定义3　BPN采用集合 $\{ (s,{{\pi}} (s))\} $ 表示 ${{M}}_i$ ，其中 ${{\pi}} (s) \in [0,1]$ 为库所 $S$ 中资源的有无，1为资源存在且满足如下规则：

1) $t_B \in T_B$ ， ${{M}}_i$ 标识处 $\forall s \in {}^ \bullet t_B:(s,t_B) \in F$ ， ${{\pi}}(s)$ 为库所获得资源的后验概率；

2)若 ${{M}}_j$ 通过变迁t的发生可直达 ${{M}}_k$ ，记为 ${{M}}_j[t > {{M}}_k$ 。

1.3 故障检测问题

定义4　给定一个 $\delta_ o \in {(L \times {\bf{R}}_0^ + )^*}$ 和一个时间 $\tau \geqslant $ $t_l(\delta _o)$ ，对应的时间序列集合与状态集合分别为

$\forall t_f \in {T}_f^i,t_f \notin {{\pi}} '$ then

$\sum {(\delta_ o,\tau )} = \{ \sigma \in {(T \times {\bf{R}}_0^ + )^*}|{{M}}_0[\sigma \rangle {{M}}\} $

$S(\delta_ o,\tau ) = \{ {{M}}_k,\varTheta _k|{{M}}_k \in {{M}})$

$\varTheta_k = \{ \bar l_{ki} \leqslant \theta_ {ki} \leqslant \bar u_{ki}\} :t_k \in {\cal{A}}({{M}}_k)\} $

定义5　给定 $\delta_ o \in {(L \times {\bf{R}}_0^ + )^*}$ 、截止时间 $\tau \in $ ${\bf{R}}_0^ + $ 和故障类型 $T_i^f$ 组成的集合 $(\delta_ o,\tau ,T_f^i)$ ，定义故障函数：

$\varGamma :{(L \times {\bf{R}}_0^ + )^*} \times {\bf{R}}_0^ + \times \{ T_f^1,T_f^2,T_f^3 \cdot \cdot \cdot T_f^i\} \to \{ N,U,F\} $

为每个集合分配一个故障状态。 $T_f^i$ 为故障类型。给定TTS序 $\sigma = (t_0,\tau_ 0)(t_1,\tau _1), \cdot \cdot \cdot, (t_h,\tau _h) \in {(T \times {\bf{R}}_0^ + )^*}$ ，记 $\log(\sigma ) = t_0,t_1,t_2, \cdot \cdot \cdot, t_h$ 为忽略触发时间 $\tau $ 的变迁序列。

1) $\forall \sigma \in \sum {(\delta _o,\tau )} ,\exists t_f \in T_f^i$ ，满足 $t_f \notin \log(\sigma )$ 则 $\varGamma (\delta_ o,\tau ,$ $T_f^i) = N$ ，故障 $t_f$ 不可能触发，任何一个TTS序列 $\sigma $ 都不包含故障变迁 $t_f$ 。

2) $\exists \sigma \in \sum {(\delta _o,\tau )} ,t_f \in T_f^i$ 满足 $t_f \in \log(\sigma )$ ，又 $\exists \sigma ' \in$ $ \sum {(\delta_ o,\tau )} ,t_f \in T_f^i$ 满足 $t_f \notin \log(\sigma )$ 。则 $\varGamma (\delta _o,\tau ,T_f^i) = {{U}}$ ，故障变迁 $t_f$ 可能触发也可能没触发，存在包含故障变迁 $t_f$ 的TTS序列 $\sigma $ ，也存在不包含故障变迁 $t_f$ 的TTS序列 $\sigma '$ 。系统发生故障概率为

$\begin{array}{l} p({\delta _o},\tau ,T_f^i) = \mathop {\max }\limits_i (p(t_f^i = 1|{t_1},{t_2}, \cdots ,{t_i})) = \\ \mathop {\max }\limits_i \Bigg(\frac{{\displaystyle p(t_f^i)\prod\limits_{a \ne i} {p({t_a}|{\rm{parents}}({t_a}))} }}{{\displaystyle\sum\limits_{k = 1}^m {[p({t_i} = {t_{ik}})\prod\limits_{a = i} {p({t_a}|{\rm{parents}}({t_a}))} ]} }}\Bigg) \end{array}$

式中： $t_1,t_2, \cdots ,t_i \in T_N$ ； $t_f^i \in T_u$ 且 $t_f^i \in T_B$ 。

3) $\forall \sigma \in \sum {(\delta_ 0,\tau )} ,\exists t_f \in T_f^i$ ，满足 $t_f \in \log(\sigma )$ ，则 $\varGamma (\delta_ o,\tau ,$ $T_f^i) = F$ ，故障 $t_f$ 肯定触发，任何一个TTS集合 $\sigma $ 都包含故障变迁 $t_f$ 。

部分可观时间Petri网故障诊断性的贝叶斯估计首先设计一个状态估计方法，遍历满足TLS序列 $\delta_ o$ 和式(6)的所有TTS序列 $\sigma $ ，分析是否存在含故障变迁 $t_f \in T_f^i$ 的TTS序列 $\sigma $ ，判断系统是否发生故障。针对POTPN方法诊断为可能故障变迁，建立对应的BPN模型，计算系统故障概率。

2 部分可观时间Petri网的贝叶斯故障检测方法 2.1 部分可观时间Petri网的状态估计

部分可观时间Petri网估计系统故障状态需要根据变迁的时间属性构建MSCG图，遍历MSCG图中满足可观变迁触发序列的路径是否包含故障变迁，估计系统故障状态。

给定一个TLS集合 $\delta_ o = (\gamma_{ i1},\tau _1)$ $(\gamma _{i2},\tau _2) ,\cdot \cdot \cdot, (\gamma_ {ih},$ $\tau_ h) \in$ $ (L \times R_0^ + )$ 和截止时间 $\tau \geqslant 0$ 对应的MSCG路径：

$ \begin{array}{c} {{\pi}} = {C_{q0}}({{{M}}_{q0}},{\varTheta _{q0}}),{t_{i1}},{\cal{L}}({t_{i1}})\\ {\varDelta ^{(1)}} \in [{\rm{max}}\{ 0,l_{i1}^{q0}\} ,{\rm{mi}}{{\rm{n}}_{j:tj}} \in {\cal{A}}({{{M}}_{qo}})\{ u_j^{q0}\} ]\\ {C_{q1}}({{{M}}_{q1}},{\Theta _{q1}}),{t_{i2}},L({t_{i2}})\\ {\varDelta ^{(2)}} \in [{\rm{max}}\{ 0,l_{i2}^{q1}\} ,{\rm{mi}}{{\rm{n}}_{j:tj}} \in {\cal{A}}({{{M}}_{q1}})\{ u_j^{q1}\} ]\\ \cdots \;\; \cdots \;\; \cdots \;\; \cdots \\ {C_{qk}}_{{\rm{ - 1}}}\;({{{M}}_{qk}}_{ - 1},{\Theta _{qk - 1}}),{t_{ik}},{\cal{L}}({t_{ik - 1}})\\ {\varDelta ^{({\rm{k}})}} \in [{\rm{max}}\{ 0,l_{ik}^{qk - 1}\} ,{\rm{mi}}{{\rm{n}}_{j:tj}} \in A({M_{qk - 1}})\{ u_j^{qk - 1}\} ]\\ {C_{qk}} \end{array} $

(5)

定理1　给定路径(5)起始于根节点 ${C_{q0}}$ 结束于 ${C_{qk}}$ ，属于可观时间序列 $\delta_ o = (\gamma_{ i1},\tau_ 1)(\gamma_{ i2},\tau_ 2) ,\cdot \cdot \cdot,$ $(\gamma_{ ih},\tau_ h) \in $ $ (L \times {\bf{R}}_0^ + )$ 对应路径集合，满足式(6)：

$\left\{ \begin{array}{l} \displaystyle\sum\limits_{{{l}} = 1}^k {{\varDelta ^{(l)}}} \leqslant \tau q \\ \tau - \displaystyle\sum\limits_{l = 1}^k {{\varDelta ^{(l)}}} < \mathop {\min \{ u_r^{qk}\} }\limits_{r:tr \in {\rm{out}}(C_{qk})} \\ {\varDelta ^{(l)}} \geqslant \max\{ 0,l_{{i_l}}^{{q_{l - 1}}}\} ,l = 1,2, \cdot \cdot \cdot ,k \\ {\varDelta ^{(l)}} \leqslant \mathop {\min \{ u_j^{{q_{l - 1}}}\} }\limits_{j:tj \in A({M_{ql}} - 1)} ,l = 1,2, \cdot \cdot \cdot ,k \\ \displaystyle\sum\limits_{l = 1}^q {{\varDelta ^{(l)}}} = \tau_ h,\forall t_{iq} \in \log (\delta_ o) \\ \end{array} \right.$

(6)

${\rm{out}}(C_q)$ 为状态集合 $C_q$ 的使能变迁集合， $\tau_ {{h}}$ 是可观测变迁 $\forall t_{iq} \in \log(\delta_ o)$ 最后触发时间。给定路径 ${{\pi}} $ 不满足式(6)则该路径不可能触发。

算法1　基于MSCG的状态估计算法。

输入　部分可观时间Petri网。

输出　故障诊断结果 $\varGamma (\delta o,\tau ,\sigma _f^i)$ 。

1) 初始化 根节点 $C_0$ ，标记 ${\cal{B}}(C) = E$ ，对应标识 ${{M}}_0$ ， $l_i^0 = l_i^{},u_i^0 = u_i^{},\varTheta 0 = \{ l_i^0 \leqslant {\theta _i} \leqslant u_i^0\} :t_i \in {\cal{A}}({M}_0)$ ；

2) while存在一个标记 ${\cal{B}}(C) = E$ 的节点do

选择一新节点 $C_k$ ，标记 ${\cal{B}}(C) = E$ ；

　for 任意 $t_i \in {\cal{A}}({{M}}_k)$ do

　　if $\max\{ 0,l_i^k\} \leqslant \min _{j:tj} \in A(M_k)\{ u_j^k\}$ then

　　　 $t_i$ 在标识 ${{M}}_k$ 触发， ${{M}}_q = {{M}}_k + C( \cdot ,t_i)$

　　for $\forall t_r \in {\cal{A}}({{M}}_q)$ do

　　if $t_r \in {\cal{A}}({{M}}_k)$ 且 ${{M}}_k - {\rm{Pre}}(:,t_i) \geqslant $ ${\rm{Pre}}(:,t_r)$

　　then

　　　 $l_r^q = l_r^k - \Delta i,u_r^q = u_r^k - \Delta i$

　　　else

　　　　 $l_r^q = l_r^{},u_r^q = u_r^{}$

建立新节点 $C_q$ 标识为 ${{M}}_q$ ，约束集合为 $\varTheta q = $ $\{ \max\{ 0,l_r^q\} \leqslant \theta_ q \leqslant u_r^q\} ,\forall t_r \in {\cal{A}}({{M}}_q)$ ， $C_q$ 和 $C_k$ 边界为：

$t_i,L(t_i),\Delta_ i \in [\max\{ 0,l_i^k\} ,$ $\min_{j : tj} \in A(M_k)\{ u_j^k\} ]$

if存在和 $C_q$ 相同节点then

$C_q$ 节点标记 ${\cal{B}}(C) = Z$

else

$C_q$ 节点标记 ${\cal{B}}(C) = E$

3) 在修正状态树(MSCT)中复制标记为 ${\cal{B}}(C) = Z$ 的节点，构建MSCG图。

4) 令 ${{II }}=\text{Ø} $ ，所有满足 ${\rm{obs}}({\pi} ) = \log(\delta_ 0)$ 且满足约束(6)和 $t_l(\sigma ) < \tau $ 的路径 $ \pi$ 放入集合 $II $ 。

5)　for $i = 1, \cdot \cdot \cdot ,r$ do

　　if $\forall {{\pi}} \;{{II}} \wedge \forall t_f \in {T}_f^i,t_f \notin {{\pi}} $ then

　　　 $\Gamma (\delta_ o,\tau ,{T}_f^i) = N$

　　else if

$\exists\;{{\pi}} \in {{II,}}{{\pi}} ' \in {{II}}$

$\exists\; t_f \in {T}_f^i,t_f \in {{\pi}} $

$\forall \;t_f \in {T}_f^i,t_f \notin {{\pi}} '$ then

　　　 $\Gamma (\delta_ o,\tau ,{ T}_f^i) = U$

　　else if $\forall {{\pi}} \in {{II,}}\forall t_f \in {T}_f^i,t_f \in {{\pi}} $

$\Gamma (\delta _o,\tau ,{T}_f^i) = F$

算法分析：算法1中步骤1)、2)、3)为构建MSCG图的过程。首先构建修正状态类树(modified state class tree，MSCT)，然后复制MSCT中标记为 ${\cal{B}}(C) =$ $ Z$ 的节点构建MSCG图。步骤4)遍历所有满足可观变迁序列 $\delta_ o$ 和式(6)的路径放在集合II中。步骤5)根据集合II中路径是否包含故障变迁判断系统故障状态。

2.2 部分可观时间Petri网的贝叶斯故障诊断

BPN模型诊断系统故障状态，首先根据前向路径中可观变迁触发状态计算故障变迁触发概率，再利用B函数根据后向路径可观变迁触发状态修正故障变迁触发概率。最后，取所有故障变迁中最大触发概率作为系统故障概率。

算法2　基于BPN模型的故障诊断算法

输入　POTPN模型及可能故障变迁 ${ T}_f^i$

输出　系统故障状态

初始化 $\forall t_f^i \in T_B,p(t_{Bm}) \in {\rm{[0,1]}}$ ；

1) for $ \forall t_f^i \in T_B $ do

　for $\forall \sigma \in \sum {(\delta_ o,\tau )} ,t_f^i \in \sigma $ do

2) While $\forall t_k \in T_o,t_k \in {{M}}_k,{{M}}_i[t_f^i\rangle , \cdots, {{M}}_j[t_j\rangle, $

　　　　 $ \cdots, {{M}}_k,\;i \leqslant j < k,\;t_j \in T_u\;\;{\rm{do}}$

复制 $t_k$ 与 $t_f^i$ 之间所有变迁 $t_j$ 、库所s和弧 $F$ ；

令 $t_j \in {T_{B,s}} \in S_B$ ，复制 $t_k \in T_o$ ，令 $t_k \in T_N$ ；

3) While $\forall_{ tk} \in T_o,t_k \in {{M}}_k,\;{{M}}_k[t_k \rangle , \cdots, {{M}}_j[t_j \rangle, $

　　　　　 $\cdots, {{M}}_i[t_j^i,k \leqslant j < i,t_j \in T_u\;{\rm{do}}$

复制 $t_k$ 与 $t_f^i$ 之间所有变迁 $t_j$ 、库所s和弧 $F$ ；

令 $t_j \in T_B,s \in S_B$ ，复制 $t_k \in T_o$ ，令 $t_k \in T_N$ ；

4) if $\exists t_k \in T_B$ 对于任意标识 ${{M}}_i$ 满足 ${{M}}_i[t_f^i \rangle $ ， $t_k$ 没有发生权 then

删除 $t_k$ 及其输入输出弧；

if $\exists s_k \in S_B$ 对于任意标识 ${{M}}_i$ 满足 ${{M}}_i[t_f^i \rangle$ ， ${{\pi}} (s_k) = 0$ then

删除 $s_k$ 及其输入输出弧；

5) for $\forall _{ti} \in T_B,\;s_i \in S_B\;{\rm{do}}$

设置 $P(t_i),\;P{\rm{(}}s_i{\rm{)}}$ ；

根据前向路径可观变迁状态计算 $p(t_f^i)$ ；

6) for $\forall t_k \in T_N,p(t_k) \to \{ 0,1\} \;{\rm{do}}$

令 $t_k \in V$ ，根据式(3)、(4)计算B函数；

7) 计算系统故障概率 $p(\delta _o,\tau ,T_f^i) = \mathop {\max}\limits_i (p(t_f^i{\rm{|}}V))$ 。

算法分析：算法2分成两部分，步骤1)、2)、3)、4)为POTPN模型基础上构建BPN模型的过程，步骤5)、6)、7)为计算系统故障概率的过程。其中步骤2)和3)根据故障变迁所在序列分别向前向后搜素第一个出现的可观变迁，并将所经过路径中不可观变迁设置为贝叶斯变迁，库所设置为贝叶斯库所，同时将可观变迁设置为常规变迁；步骤4)根据贝叶斯变迁与故障变迁是否存在竞争关系，删除不可能触发贝叶斯变迁及对应的贝叶斯库所，降低模型复杂度；步骤5)计算故障变迁触发的条件概率；步骤6)根据后向路径可观变迁触发状态修正故障变迁触发概率，然后取所有故障变迁中最大触发概率作为系统故障(步骤7))。

2.3 算法复杂分析

本文故障诊断方法包含算法1和算法2两部分，下面分析算法1的计算复杂度。

算法1首先构建MSCG图，根据可观变迁触发序列遍历所有路径，分析是否存在包含故障变迁的路径，判断系统故障状态。其中步骤2)和3)根据变迁 $t_i \in {\cal{A}}({{M}}_q)$ 在标识 ${{M}}_k$ 使能情况，修正时间约束，建立新节点 $C_q$ ，最坏情况变迁 $t_i$ 每触发一次，所有变迁都使能，对应计算复杂度为 ${\cal{O}}({\left| T \right|^2})$ 。给定截止时间 $\tau $ ，满足 $\delta_ 0$ 的路径个数有限，定义 $l_{{\rm{max}}}$ 为MSCG中满足 $\delta_ 0$ 的最长路径，满足 $\delta_ 0$ 的 $\sigma $ 数最多为 ${\left| T \right|^{l_{\rm{max}}}}$ ，即步骤2)、3)执行次数最多为 ${\left| T \right|^{l_{\rm{max}}}}$ ，对应计算复杂度为 ${\cal{O}}({\left| T \right|^{l_{\rm{max}}}})$ 。步骤4)、5)判断所有满足 ${\rm{obs}}({{\pi}} ) = \varepsilon $ 路径及 ${\rm{obs}}({{\pi}} ) = \delta_ 0$ 的路径是否满足式(6)并将路径放入相应集合，对应计算复杂度为 ${\cal{O}}(l_{{\rm{max}}}{\left| T \right|^{l_{\rm{max}}}})$ 。综上，算法1的计算复杂度为 ${\cal{O}}(l_{{\rm{max}}}{\left| T \right|^{l_{\rm{max}}}})$ 。

算法2首先把所有故障变迁放入 $T_B$ ，令故障变迁个数为N，步骤2)、3)以每个故障变迁为开始，分别向前向后遍历第一个可观变迁，最差时每个都要遍历 $\dfrac{{l_{{\rm{max}}}}}{2}$ ， $l_{\rm{max}} $ 为MSCG中满足 $\delta_0 $ 的最长路径，N个故障变迁则为 $\dfrac{{l_{{\rm{max}}}}}{2} N$ ，即此时复杂度为 ${\left( {\dfrac{{l_{{\rm{max}}}}}{2} N} \right)^2}$ ，第4)步去掉无关库所，即遍历所有 ${\pi} (s_k) \in$ $ S_B$ ，令数量为S，则复杂度为S，第5)步设置贝叶斯库所和变迁(N+S)触发概率，步骤6)、7)计算故障概率，复杂度为贝叶斯变迁数量N，则算法2的复杂度为 ${N^2}\left( {N + S} \right){\left( {\dfrac{{l_{{\rm{max}}}}}{2} N} \right)^2}$ 。

3 仿真分析与验证 3.1 部分可观时间Petri网的贝叶斯故障诊断 3.1.1 航天推进系统POTPN模型构建

双组元挤压推进系统工作原理为高压氦气通入推进剂储箱挤压推进剂进入下游燃烧室混合点火产生推力^[16-17]。图1以推进系统工作过程中关键节点为库所，关键动作为变迁建立Petri网模型，模拟推进系统工作过程，各库所、变迁含义如表1、2所示。

双组元推进系统工作过程为：液路管路预冷完成后地面控制站发出动作命令，气路电爆阀起爆，减压阀和单向阀受气体压力影响开启，储箱开始增压，液态燃料管路阀门和液态氧化剂管路阀门相继打开，推进剂在燃烧室混合点火，发动机工作。

3.1.2 航天推进系统故障状态估计

双组元推进系统中加注阀故障会导致推进系统故障，影响系统正常工作。本文Petri网模型中，无故障情况下燃料加注阀（变迁t₅）会在燃料储箱增压后3.1~5.2 ms内打开，若电磁阀未打开(变迁t₈)或燃料加注阀泄露(变迁t₉)，则系统发生故障。

第1次实验观测序列为 $(t_0,0)(t_3,7.6)(t_4,11.7)$ $(t_{13},24.9)$ ，其中 $(t_0,0) (t_3,7.6)$ $(t_4,11.7)$ 触发后，计算所有满足时间约束的TTS集合σ共3条，正常变迁σ为3条，故障变迁σ为0条，最大同步距离 $s_d\left( {{t_3},{t_4}} \right)$ 为2，表明系统目前运行正常。可观测变迁 触发后，计算所有满足时间约束的TTS集合共9条，其中正常变迁为3条，故障变迁 为6条，最大同步距离为2，电磁阀未打开故障和燃料加注阀泄露故障都有可能触发。

第2次实验，观测序列为 $(t_0,0)(t_4,8.1)(t_3,8.7)$ $(t_{13},29.1)$ ，其中 $(t_0,0)(t_4,8.1)$ $(t_3,8.7)$ 触发后，计算所有满足时间约束的TTS集合 $\sigma $ 共2条，正常变迁 $\sigma $ 为2条，故障变迁 $\sigma $ 为0条，最大同步距离 $s_d(t_3,t_4)$ 为2，表明系统目前运行正常。可观测变迁 $(t_{13},29.1)$ 触发后，计算所有满足时间约束的TTS集合 $\sigma $ 共4条，其中正常变迁 $\sigma $ 为2条，故障变迁 $\sigma $ 为2条，最大同步距离 $s_d(t_3,t_4)$ 为2，燃料加注阀泄露故障可能触发。

3.2 基于部分可观时间Petri网的航天推进系统贝叶斯故障诊断 3.2.1 航天推进系统BPN模型构建

根据POTPN故障诊断方法，2次实验诊断出系统故障状态都为可能故障，其中第1次实验可能故障变迁为t₈（电磁阀未打开故障）和t₉（燃料加注阀故障），第2次实验可能故障变迁为t₉（燃料加注阀故障）。利用BPN故障诊断方法判断系统故障状态，首先需要根据算法2分别构建故障变迁t₈和故障变迁t₉对应的BPN模型，然后计算2个故障变迁各自触发概率以及系统故障概率。

在推进系统的部分可观时间Petri网模型中，给定故障变迁，向前向后搜索其所有路径中的变迁和库所，直到出现可观测变迁为止。对于故障变迁 $t_9$ ，其所在路径包括 $t_3t_9t_{11},t_3t_9t_{10},t_3t_9t_7t_6t_{13}t_{12}$ ，其中路径 $t_3t_9t_{11}$ 中的可观变迁为 $t_3,t_{11}$ ，不可观变迁为 $t_9$ ；路径 $t_3t_9t_{10}$ 中可观变迁为 $t_3,t_{10}$ ，不可观变迁为 $t_9$ ；路径 $t_3t_9t_7t_6t_{13}t_{12}$ 中可观变迁为 $t_3,t_{13},t_{12}$ ，不可观变迁为 $t_9t_7t_6t$ 。令 $t_9,t_7,t_6$ 为贝叶斯变迁，库所 $s_5,s_6,s_7,s_8,s_9,$ $s_{10},s_{11},s_{12}$ 为贝叶斯库所，建立 $t_9$ 燃料加注阀泄露故障对应的BPN模型诊断子网如图2所示。

3.2.2 基于BPN的故障诊断结果

为了在BPN模型中使用贝叶斯推理计算系统故障概率，需获得贝叶斯变迁和库所的条件概率表及可观变迁触发状态，这些信息体现了故障过程的不确定性。本文基于POTPN的航天推进系统故障诊断中，2次仿真实验结果为可能故障，需要利用BPN故障诊断算法确定系统故障概率。

表3是BPN算法诊断结果，第一次试验中，待诊断故障变迁为t₉和t₈，监控中心收到如下可观变迁触发状态：可观变迁t₁₀、t₁₁、t₁₂、t₁₃触发。首先计算故障变迁t₉ 触发概率，由可观变迁t₃触发，计算故障变迁t₉触发的条件概率为0.018 7， ${{\pi}}(s_5) = 0.982$ ；根据可观变迁 $t_{12},t_{13}$ 触发状态和式(3)计算B函数，修正故障变迁t₉的触发概率为0.018 9，对应路径中贝叶斯库所获得资源概率为 ${{\pi}}(s_9) = 1,{{\pi}}(s_{12}) = 1$ , ${{\pi}}(s_8) = 0.987,{{\pi}}(s_7) = 0.996$ ；根据可观变迁t₁₀触发状态和式(3)计算B函数，修正故障变迁t₉的触发概率为0.074 9，对应路径中贝叶斯库所获得资源概率为 ${{\pi}}(s_{10}) = 1$ ；根据可观变迁t₁₁触发状态和式(3)计算B函数，修正故障变迁t₉的触发概率为0.963 0，对应路径贝叶斯库所获得资源概率为 ${{\pi}}(s_{11}) = 1$ 。同理故障变迁t₈触发概率为，则 ${\rm{max}}(p({t_8} = 1 | {t_3} = 1,{t_{12}} = 1,{t_{13}} = 1),$ $p({t_9} = 1 |{t_3} = 1, $ ${t_{10}} = 1,{t_{11}} = 1,t_{12} = 1,t_{13} = 1)) = 0.963\;0$ ，系统故障概率为0.963 0。

第2次试验中，待诊断故障变迁为t₉，监控中心收到如下可观变迁触发状态：t₁₀、t₁₁、t₁₂未触发，t₁₃触发。由可观变迁t₃触发，计算故障变迁t₉触发的条件概率为0.018 7， ${\pi} (s_5) = 0.982$ ；根据可观变迁触发 $t_{12},t_{13}$ 状态和式(3)计算B函数，修正故障变迁t₉的触发概率为0.018 7，对应路径贝叶斯库所获得资源概率 ${{\pi}} ({s_9}){\rm{ = }}1,$ ${{\pi}} (s_{12}){\rm{ = }}0.040\;2,$ ${{\pi}} ({s_8}){\rm{ = }}0.954,$ ${{\pi}} (s_7){\rm{ = }}$ $0.969$ ；根据可观变迁t₁₀触发状态和式(3)计算B函数，修正故障变迁t₉的触发概率为0.003 6，对应路径贝叶斯库所获得资源概率为 ${{\pi}} ({s_{10}}) = 0.082$ ；根据可观变迁t₁₁触发状态和式(3)计算B函数，修正故障变迁t₉的触发概率为0.008，对应路径中贝叶斯库所获得资源概率为 ${{\pi}} ({s_{11}}) = 0.035$ 。系统只存在一个可能故障变迁， ${\rm{max}}(p(t_9{\rm{ = 1|}}t_3{\rm{ = 1,}}t_{10}{\rm{ = }}0{\rm{,}}t_{11}{\rm{ = }}0{\rm{,}}$ $t_{12}{\rm{ = }}0{\rm{,}}t_{13}{\rm{ = 1}})){\rm{ = }}0.008$ ，系统故障概率为0.008。

4 结论

1)本文基于变迁时间约束，提出一种构建MSCG图的故障诊断方法，根据可观序列集合和MSCG进行系统故障诊断。

2)建立故障变迁对应的BPN模型，根据可观变迁触发状态计算不可观变迁触发概率和系统故障触发概率。

3)建立双组元推进系统的POTPN网络故障诊断模型，根据系统工作原理和各部分触发关系修正变迁时间区间，构建对应的BPN模型，通过本研究提出的算法诊断系统的故障状态。本研究提出的故障诊断算法能够满足简单系统的故障诊断需要。

之后的研究，将考虑诊断算法的诊断效率问题。当系统关键的库所和变迁的个数过多时，MSCG图的复杂性会快速上升，对于复杂性过高的系统，该方法需进一步改善。