CN104850750A - 一种核电站反应堆保护系统可靠性分析方法 - Google Patents

Abstract

本发明提供一种核电站反应堆保护系统可靠性分析方法，包括根据核电站反应堆保护系统的结构及功能联系确定其故障树模型；求解引发顶事件的最小割集；获取核电站反应堆保护系统的历史故障数据；计算核电站反应堆保护系统的寿命统计量；建立核电站反应堆保护系统寿命分布概率密度函数的信息熵模型，求解核电站反应堆保护系统最优的寿命分布概率密度函数、核电站反应堆保护系统的失效概率、核电站反应堆保护系统可靠度、核电站反应堆保护系统动态失效率。本发明依据少量可靠性试验数据对系统整体的寿命分布与动态失效率做出预测，与基于大样本的蒙特卡洛模拟结果相符，为小失效概率情况下核电站反应堆保护系统寿命预测和动态失效率评估提供技术方法。

Description

一种核电站反应堆保护系统可靠性分析方法

技术领域

本发明属于系统可靠性分析与设计技术领域，具体涉及一种核电站反应堆保护系统可靠性分析方法。

背景技术

系统可靠性表示系统在规定的条件下和规定的时间内完成规定功能的能力。从整个全生命周期的视角看，系统能否能够完成预期功能有多个衡量指标：对于可修系统和设备，衡量指标包括可靠度、平均故障间隔时间(Mean Time Between Failure，MTBF)、平均修复时间(MeanTime To Repair，MTTR)、可用度、有效寿命等；对于不可修系统或产品，包括可靠度、可靠寿命、故障率、平均寿命(Mean Time To Failure，MTTF)等技术指标。

产品设计结束后需要经过严格的材料物理试验与应力筛选、生产和制造过程的工艺控制，以及严格的质量检测等环节，这时产品所具有的可靠性被称为“固有可靠性”。可靠性试验是产品研制单位和使用单位了解产品可靠性、获得可靠性数据的基本途径。由于产品寿命试验具有破坏性，设计和生产单位一般通过少量的试验数据预测产品总体的寿命可靠性水平和各类可靠性指标，这需要一种依据小子样能够对系统可靠性水平和动态失效率做出正确预计的高精度技术方法。

核电站反应堆保护系统的失效事件具有影响后果严重、影响时间长、影响范围大等高风险的特点，所以系统各个工作部件和结构的可靠性标准要求很高，失效概率一般在十万次一遇或百万次一遇的水平。应用传统可靠性估计方法，理论上需要每百万次的可靠性试验才能获得一个失效数据样本，必然产生巨大的人力、物力和财力等成本消耗。如果能够提供一种基于几百次或千次小子样可靠性试验，能够获得10^-5～10^-6级别系统失效概率的精确估计技术方法，无疑对解决复杂机电系统小失效概率的评价问题具有重要的工程意义。

发明内容

针对现有技术存在的问题，本发明提供一种核电站反应堆保护系统可靠性分析方法。

本发明的技术方案是：

一种核电站反应堆保护系统可靠性分析方法，包括以下步骤：

步骤1、根据核电站反应堆保护系统的结构及功能联系确定其故障树模型；

步骤2、用下行法求解引发顶事件即核电站紧急停堆失效发生的最小割集，即底事件同时发生则导致顶事件发生的组合；

步骤3、获取核电站反应堆保护系统的历史故障数据即各个底事件发生故障的时间统计数据样本；

步骤4、通过最小割集计算核电站反应堆保护系统的寿命统计量，即核电站反应堆保护系统发生故障时已经运行的时间；

步骤5、建立核电站反应堆保护系统寿命分布概率密度函数的信息熵模型，求解核电站反应堆保护系统最优的寿命分布概率密度函数；

步骤6、利用核电站反应堆保护系统寿命分布概率密度函数求解核电站反应堆保护系统的失效概率，即核电站反应堆保护系统在t时刻之前发生故障的概率；

步骤7、利用核电站反应堆保护系统寿命分布概率密度函数求解核电站反应堆保护系统可靠度，即核电站反应堆保护系统在时间t之后仍然正常工作的概率；

步骤8、利用核电站反应堆保护系统寿命分布概率密度函数求解核电站反应堆保护系统动态失效率，即核电站反应堆保护系统工作到时刻t时尚未失效、在时刻t以后的单位时间内发生失效的概率。

所述步骤1包括以下步骤：

步骤1.1、确定故障树的顶事件及导致顶事件发生的二级事件；

故障树的顶事件为核电站紧急停堆失效；

导致顶事件的二级事件包括稳压器压力低信号失效、停堆断路器拒开、至少三束控制棒卡住；

任一二级事件发生均导致顶事件发生；

步骤1.2、确定故障树中导致二级事件发生的三级事件；

导致稳压器压力低信号失效事件发生的事件包括稳压器压力传感器失效、三个稳压器压力传感器阈值继电器定值错误；

导致停堆断路器拒开事件发生的事件为两个停堆断路器共因失效；

至少三束控制棒卡住事件视为底事件即导致核电站紧急停堆失效的不能再分的事件；

步骤1.3、确定导致三级事件发生的可能事件即四级事件；

导致稳压器压力传感器失效的事件包括相关的稳压器压力传感器失效；其中任意两个事件同时发生则会导致稳压器压力传感器失效事件的发生；

三个稳压器压力传感器阈值继电器定值错误、两个停堆断路器共因失效均为底事件；

步骤1.4、确定导致四级事件发生的可能事件即五级事件；

步骤1.5、确定导致五级事件发生的可能事件即六级事件，直至该事件不能再分为止，引发顶事件即核电站紧急停堆失效事件发生的所有可能事件搜寻完毕，得到底事件。

步骤5所述核电站反应堆保护系统寿命分布概率密度函数按如下步骤建立：

步骤5.1、引入最大熵估计方法，建立核电站反应堆保护系统的信息熵模型；

步骤5.2、确定核电站反应堆保护系统信息熵优化模型的约束条件，包括核电站反应堆保护系统寿命统计量的分数矩的建立、核电站反应堆保护系统寿命分布概率密度函数估计的积分值为1；

步骤5.3、引入拉格朗日方程求解核电站反应堆保护系统最大熵约束优化问题，使拉格朗日方程对核电站反应堆保护系统寿命分布概率密度估计量求偏导数，令其值等于0，得到核电站反应堆保护系统寿命分布概率密度函数的解析式估计；

步骤5.4、引入K-L距离方法，建立求解核电站反应堆保护系统寿命分布概率密度函数参数的无约束优化模型；

步骤5.5、利用历史故障数据求解核电站反应堆保护系统寿命分布概率密度函数的拉格朗日乘子λ与分数矩指数α；

步骤5.6、将拉格朗日乘子λ与分数矩指数α代入步骤5-3中的核电站反应堆保护系统的概率密度函数的解析式估计，得到核电站反应堆保护系统寿命分布概率密度函数。

有益效果：

本发明突破了传统方法对大量观测样本的依赖，依据少量可靠性试验数据对系统整体的寿命分布与动态失效率做出预测，计算结果与基于大样本的蒙特卡洛模拟结果相符，为小失效概率情况下核电站反应堆保护系统的寿命预测和动态失效率评估提供技术方法。

附图说明

图1是本发明具体实施方式的核电站反应堆保护系统可靠性分析方法流程图；

图2是本发明具体实施方式的核电站反应堆保护系统故障树模型；

图3是本发明具体实施方式的简化后的故障树模型；

图4是本发明具体实施方式的核电站反应堆保护系统寿命分布概率密度函数建立流程图；

图5是本发明具体实施方式的核电站反应堆保护系统寿命分布概率密度函数；

图6是本发明具体实施方式的核电站反应堆保护系统失效概率预测曲线；

图7是本发明具体实施方式的核电站反应堆保护系统可靠度预测曲线；

图8是本发明具体实施方式的核电站反应堆保护系统动态失效率预测曲线。

具体实施方式

下面结合附图对本发明的具体实施方式做详细说明。

以某核电站反应堆保护系统的可靠性分析为例，详细说明本发明方法的实施过程，该核电站反应堆保护系统可靠性分析方法，如图1所示，包括以下步骤：

步骤1、根据核电站反应堆保护系统的结构及功能联系确定其故障树模型；

步骤1.1、确定故障树的顶事件及导致顶事件发生的二级事件；

故障树的顶事件为核电站紧急停堆失效(RCPS000)；

导致顶事件的二级事件包括稳压器压力低信号失效(RCPS001)、停堆断路器拒开(RCPS002)、至少三束控制棒卡住(RCPS003)；

任一二级事件发生均导致顶事件(RCPS000)发生；

步骤1.2、确定故障树中导致二级事件发生的三级事件；

导致稳压器压力低信号失效事件(RCPS001)发生的事件包括稳压器压力传感器失效(RCPS004)、三个稳压器压力传感器阈值继电器定值错误(PCF005-013)；

导致停堆断路器拒开事件(RCPS002)发生的事件为两个停堆断路器共因失效(RPA300JA-RO)；

至少三束控制棒卡住事件(RCPS003)视为底事件即导致核电站紧急停堆失效的不能再分的事件；

步骤1.3、确定导致三级事件发生的可能事件即四级事件；

导致稳压器压力传感器失效(RCPS004)的事件包括相关的稳压器压力传感器失效，具体是稳压器压力传感器RCP005MP失效(RCPS005)、稳压器压力传感器RCP006MP失效(RCPS006)和稳压器压力传感器RCP007MP失效(RCPS007)；其中任意两个事件同时发生则会导致稳压器压力传感器失效(RCPS004)事件的发生；

三个稳压器压力传感器阈值继电器定值错误(PCF005-013)、两个停堆断路器共因失效(RPA300JA-RO)均为底事件；

步骤1.4、确定导致四级事件发生的可能事件即五级事件；

导致稳压器压力传感器RCP005MP失效(RCPS005)事件发生的事件有稳压器压力传感器RCP005MP需求失效(RCPS005-MP)、稳压器压力传感器RCP005MP阈值继电器失效(RCPS005-RC)，其中任意事件发生均会导致稳压器压力传感器RCP005MP失效(RCPS005)事件发生。

导致稳压器压力传感器RCP006MP失效(RCPS006)事件发生的事件有稳压器压力传感器RCP006MP需求失效(RCPS006-MP)和稳压器压力传感器RCP006MP阈值继电器失效(RCPS006-RC)，其中任意事件发生均会导致稳压器压力传感器RCP006MP失效(RCPS006)事件发生。

导致稳压器压力传感器RCP007MP失效(RCPS007)事件发生的事件有稳压器压力传感器RCP007MP需求失效(RCPS007-MP)和稳压器压力传感器RCP007MP阈值继电器失效(RCPS007-RC)，其中任意事件发生均会导致稳压器压力传感器RCP007MP失效(RCPS007)事件发生。

步骤1.5、确定导致五级事件发生的可能事件即六级事件，直至该事件不能再分为止，引发顶事件即核电站紧急停堆失效事件发生的所有可能事件搜寻完毕，得到底事件。

稳压器压力传感器RCP005MP需求失效(RCPS005-MP)、稳压器压力传感器RCP005MP阈值继电器失效(RCPS005-RC)、稳压器压力传感器RCP006MP需求失效(RCPS006-MP)、稳压器压力传感器RCP006MP阈值继电器失效(RCPS006-RC)、稳压器压力传感器RCP007MP需求失效(RCPS007-MP)和稳压器压力传感器RCP007MP阈值继电器失效(RCPS007-RC)均为底事件。

至此引发顶事件即核电站紧急停堆失效(RCPS000)事件发生的所有可能事件搜寻完毕。可根据其父子关系建立该系统的故障树模型，如图2所示。

步骤2、用下行法求解引发顶事件即核电站紧急停堆失效(RCPS000)发生的最小割集，即底事件同时发生则导致顶事件发生的组合；

为简化计算，事件核电站紧急停堆失效(RCPS000)用T表示，事件稳压器压力低信号失效(RCPS001)用G1表示，事件停堆断路器拒开(RCPS002)用G2表示，事件至少三束控制棒卡住(RCPS003)用A表示，事件稳压器压力传感器失效(RCPS004)用G3表示，事件三个稳压器压力传感器阈值继电器定值错误(PCF005-013)用B表示，两个事件停堆断路器失效(RPA300JA-RO)分别用C和D表示，事件稳压器压力传感器RCP005MP失效(RCPS005)用G4表示，事件稳压器压力传感器RCP006MP失效(RCPS006)用G5表示，事件稳压器压力传感器RCP007MP失效(RCPS007)用G6表示，事件稳压器压力传感器RCP005MP需求失效(RCPS005-MP)用E表示，事件稳压器压力传感器RCP005MP阈值继电器失效(RCPS005-RC)用F表示，事件稳压器压力传感器RCP006MP需求失效(RCPS006-MP)用G表示、事件稳压器压力传感器RCP006MP阈值继电器失效(RCPS006-RC)用H表示、事件稳压器压力传感器RCP007MP需求失效(RCPS007-MP)用I表示，事件稳压器压力传感器RCP007MP阈值继电器失效(RCPS007-RC)用J表示，简化后的故障树模型如图3所示。

求解故障树的最小割集的过程如下：

第一步：顶事件T下面的门为或门(任一事件发生，顶事件即发生)，因此将它的输入G1、G2、A排成列(置换T)；

第二步：基本事件A不再分解，G1事件下为或门，将其输入G3、B排成一列置换G1；G2事件下为与门(所有事件都发生，顶事件才发生)，将其输入C、D排成一行，置换G2；

第三步：基本事件B、C、D不再分解，G3事件下为2/3表决门(子事件中有两个或两个以上发生，顶事件发生)，将其输入G4、G5、G6两两组合成行后，排成列置换G3；

第四步：G4事件下为或门，将其输入E、F排成一列置换G4,；G5事件下为或门，将其输入G、H排成一列置换G5；G6事件下为或门，将其输入I、J排成一列置换G6；

第五步：得到一列全部由基本事件表示的15个割集：{A}，{CD}，{B}，{EG}，{EH}，{FG}，{FH}，{EI}，{EJ}，{FI}，{FJ}，{GI}，{GJ}，{HI}，{HJ}，见表1。

表1核电站反应堆保护系统故障树下行法展开

综上，可得到致使核电站紧急停堆失效(RCPS000)发生的最小的底事件组合。

步骤3、获取核电站反应堆保护系统的历史故障数据即各个底事件(A，B，C，D，E，F，G，H，I，J)发生故障的时间统计数据样本：

$\{t_A^{\left(1\right)},t_A^{\left(2\right)},...,t_A^{\left(N\right)}\},\{t_B^{\left(1\right)},t_B^{\left(2\right)},...,t_B^{\left(N\right)}\},\{t_C^{\left(1\right)},t_C^{\left(2\right)},...,t_C^{\left(N\right)}\},\{t_D^{\left(1\right)},t_D^{\left(2\right)},...,t_D^{\left(N\right)}\},\{t_E^{\left(1\right)},t_E^{\left(2\right)},...,t_E^{\left(N\right)}\},$

$\{t_F^{\left(1\right)},t_F^{\left(2\right)},...,t_F^{\left(N\right)}\},\{t_G^{\left(1\right)},t_G^{\left(2\right)},...,t_G^{\left(N\right)}\},\{t_H^{\left(1\right)},t_H^{\left(2\right)},...,t_H^{\left(N\right)}\},\{t_I^{\left(1\right)},t_I^{\left(2\right)},...,t_I^{\left(N\right)}\},\{t_J^{\left(1\right)},t_J^{\left(2\right)},...,t_J^{\left(N\right)}\}.$

其中，t表示发生故障时已经正常运行的时间，上标表示不同的样本，N为样本数量，下标表示不同的事件。

步骤4、通过最小割集计算核电站反应堆保护系统的寿命统计量，即核电站反应堆保护系统发生故障时已经运行的时间

第i个时间统计数据样本为最小割集中最少的工作时间，含有两个或两个以上底事件元素的工作时间为其中最大的工作时间。具体的为：从第i个时间统计数据样本中，找出至少三束控制棒卡住(RCPS003)即事件A发生的时间，三个稳压器压力传感器阈值继电器定值错误(PCF005-013)即事件B发生的时间，停堆断路器失效(RPA300JA-RO)即事件C和D发生的时间的较大者，稳压器压力传感器RCP005MP需求失效(RCPS005-MP)和稳压器压力传感器RCP006MP需求失效(RCPS006-MP)即事件E和事件G发生时间的较大者，稳压器压力传感器RCP005MP需求失效(RCPS005-MP)和稳压器压力传感器RCP006MP阈值继电器失效(RCPS006-RC)即事件E和H发生时间的较大者，稳压器压力传感器RCP005MP阈值继电器失效(RCPS005-RC)和稳压器压力传感器RCP006MP需求失效(RCPS006-MP)即时间F和G发生时间的较大者，稳压器压力传感器RCP005MP阈值继电器失效(RCPS005-RC)和稳压器压力传感器RCP006MP阈值继电器失效(RCPS006-RC)即事件F和H发生时间的较大者，稳压器压力传感器RCP005MP需求失效(RCPS005-MP)和稳压器压力传感器RCP007MP需求失效(RCPS007-MP)即事件E和I发生时间的较大者，稳压器压力传感器RCP005MP需求失效(RCPS005-MP)和稳压器压力传感器RCP007MP阈值继电器失效(RCPS007-RC)即事件E和J发生时间的较大者，稳压器压力传感器RCP005MP阈值继电器失效(RCPS005-RC)和稳压器压力传感器RCP007MP需求失效(RCPS007-MP)即事件F和I发生时间的较大者，稳压器压力传感器RCP005MP阈值继电器失效(RCPS005-RC)和稳压器压力传感器RCP007MP阈值继电器失效(RCPS007-RC)即事件F和J发生时间的较大者，稳压器压力传感器RCP006MP需求失效(RCPS006-MP)和稳压器压力传感器RCP007MP需求失效(RCPS007-MP)即事件G和I发生时间的较大者，稳压器压力传感器RCP006MP需求失效(RCPS006-MP)和稳压器压力传感器RCP007MP阈值继电器失效(RCPS007-RC)即事件G和J发生时间的较大者，稳压器压力传感器RCP006MP阈值继电器失效(RCPS006-RC)和稳压器压力传感器RCP007MP需求失效(RCPS007-MP)即事件H和I发生时间的较大者，稳压器压力传感器RCP006MP阈值继电器失效(RCPS006-RC)和稳压器压力传感器RCP007MP阈值继电器失效(RCPS007-RC)即事件H和J发生时间的较大者，从以上时间中取最小的时间即为第i个时间统计数据样本的工作时间。

其具体可以通过下式计算：

$t_s^{\left(i\right)}=\min \left\{\begin{array}{c}{t}_A^{\left(i\right)},t_B^{\left(i\right)},\max \{t_C^{\left(i\right)},t_D^{\left(i\right)}\},\max \{t_E^{\left(i\right)},t_G^{\left(i\right)}\},\max \{t_E^{\left(i\right)},t_H^{\left(i\right)}\},\max \{t_F^{\left(i\right)},t_G^{\left(i\right)}\},\\ \max \{t_F^{\left(i\right)},t_H^{\left(i\right)}\},\max \{t_E^{\left(i\right)},t_I^{\left(i\right)}\},\max \{t_E^{\left(i\right)},t_J^{\left(i\right)}\},\max \{t_F^{\left(i\right)},t_I^{\left(i\right)}\},\max \{t_F^{\left(i\right)},t_J^{\left(i\right)}\},\\ \max \{t_G^{\left(i\right)},t_I^{\left(i\right)}\},\max \{t_G^{\left(i\right)},t_J^{\left(i\right)}\},\max \{t_H^{\left(i\right)},t_I^{\left(i\right)}\},\max \{t_H^{\left(i\right)},t_J^{\left(i\right)}\},\end{array}\right\};$

由此依次可得出系统的寿命统计量

步骤5、建立核电站反应堆保护系统寿命分布概率密度函数的信息熵模型，求解核电站反应堆保护系统最优的寿命分布概率密度函数；

如图4所示，核电站反应堆保护系统寿命分布概率密度函数按如下步骤建立：

步骤5.1、引入最大熵估计方法，建立核电站反应堆保护系统的信息熵模型；

为了确定寿命分布中的分数矩指数α和拉格朗日乘子λ，则需要引入最大熵估计方法。

已知寿命分布概率密度函数f_T(t)，其信息熵模型定义为：

H[f]＝-∫_Tf_T(t)log[f_T(t)]dt

步骤5.2、确定核电站反应堆保护系统信息熵优化模型的约束条件，包括核电站反应堆保护系统寿命统计量的分数矩的建立、核电站反应堆保护系统寿命分布概率密度函数估计的积分值为1；

随机变量T为部件(或系统)寿命随机变量，其α阶分数矩定义为：

$M_T^{\mathrm{\α}}=E\left[T^{\mathrm{\α}}\right]={\∫}_T{t}^{\mathrm{\α}}{f}_T\left(t\right)\mathrm{dt}$

其中α为任意实数。

需要说明的是：寿命随机变量分数矩存在的充分必要条件是分数矩积分收敛，其等价于存在k阶整数矩，当且仅当|α|≤k时存在。

为了计算寿命分布概率密度函数fT(t)的估计量则需要应用最大熵理论。通过引入寿命统计量的各阶分数矩，信息熵优化模型的约束条件为： ${\∫}_T{t}^{{\mathrm{\α}}_j}{\hat{f}}_T\left(t\right)\mathrm{dt}=M_T^{{\mathrm{\α}}_j},(j=\mathrm{1,2},...,m)$ 和 ${\∫}_T{\hat{f}}_T\left(t\right)\mathrm{dt}=1,$ 其中m为分数矩次数，经多次实践应用验证，取三阶分数矩约束(即m＝3)即可达到获得满意计算精度；αj为对应阶次的分数矩指数；为核电站反应堆保护系统寿命即故障前工作时间的分数矩估计，N为系统的时间统计数据样本数量。此时，只需要使信息熵估计最大化即可。

其信息熵优化模型可以表示为：

$\left\{\begin{array}{c}\mathrm{Find}:{\hat{f}}_T\left(t\right)\\ \mathrm{Maximize}:H\left[\hat{f}\right]=-{\∫}_T{\hat{f}}_T\left(t\right)\log \left[{\hat{f}}_T\left(t\right)\right]\mathrm{dt}\\ \mathrm{Subject\; to}:{\∫}_T{t}^{{\mathrm{\α}}_j}{\hat{f}}_T\left(t\right)\mathrm{dt}=M_T^{{\mathrm{\α}}_j},(j=\mathrm{1,2},...,m)\end{array}\right.$

步骤5.3、引入拉格朗日方程求解核电站反应堆保护系统最大熵约束优化问题，使拉格朗日方程对核电站反应堆保护系统寿命分布概率密度估计量求偏导数，令其值等于0，得到核电站反应堆保护系统寿命分布概率密度函数的解析式估计；

引入拉格朗日方程求解该最大熵约束优化问题：

$\begin{array}{c}L[\mathrm{\λ},\mathrm{\α};{\hat{f}}_T\left(t\right)]=-{\∫}_T{\hat{f}}_T\left(t\right)\log \left[{\hat{f}}_T\left(t\right)\right]\mathrm{dt}-({\mathrm{\λ}}_0-1)[{\∫}_T{\hat{f}}_T\left(t\right)\mathrm{dt}-1]\\ -\underset{j=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\λ}}_j[{\∫}_T{t}^{{\mathrm{\α}}_j}\left(t\right)\mathrm{dt}-M_T^{{\mathrm{\α}}_j}]\end{array}$

其中λ＝[λ₀,λ₁,…,λ_m]^T为拉格朗日乘子，α＝[α₀,α₁,…,α_m]^T为分数矩指数。

为得到最大熵估计，只需令拉格朗日方程对概率密度估计量求偏导数令其等于0即可。即：令 $\frac{\mathrm{\δL}[\mathrm{\λ},\mathrm{\α};{\hat{f}}_T\left(t\right)]}{\mathrm{\δ}{\hat{f}}_T\left(t\right)}=0$

得到未知寿命概率密度函数的分数矩信息熵估计表达式：

${\hat{f}}_T\left(t\right)=\exp (-\underset{j=0}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\λ}}_j{t}^{{\mathrm{\α}}_j}),({\mathrm{\α}}_0=0)$

考虑到求得拉格朗日乘子λ₀的表达式为：

${\mathrm{\λ}}_0=\log \left[{\∫}_T\exp (-\underset{j=0}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\λ}}_j{t}^{{\mathrm{\α}}_j})\mathrm{dt}\right]$

为了简化信息熵优化模型，引入K-L距离求解寿命分布概率密度函数的拉格朗日乘子λ与分数矩指数α。

步骤5.4、引入K-L距离方法，建立求解核电站反应堆保护系统寿命分布概率密度函数参数的无约束优化模型；

K-L距离定义为真实熵和估计熵之间的差值，其越小表明估计熵越接近真实熵，就越准确。其公式为：

$\begin{array}{c}K[f,\hat{f}]={\∫}_T{f}_T\left(t\right)\log [f_T\left(t\right)/{\hat{f}}_T\left(t\right)]\mathrm{dt}\\ ={\∫}_T{f}_T\left(t\right)\log \left[f_T\left(t\right)\right]\mathrm{dt}-{\∫}_T\log \left[{\hat{f}}_T\left(t\right)\right]f_T\left(t\right)\mathrm{dt}\end{array}$

当给定寿命分布概率密度函数估计量的表达式K-L距离可以进一步表示为：

$K[f,\hat{f}]=-H\left[f\right]+{\mathrm{\λ}}_0+\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\λ}}_i{M}_T^{{\mathrm{\α}}_i}$

由于H[f]为寿命分布的信息熵的理论值通常为实常数。因此K-L距离的变化部分可以表示为：

$I(\mathrm{\λ},\mathrm{\α})=\log \left[{\∫}_T\exp (-\underset{j=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\λ}}_j{t}^{{\mathrm{\α}}_j})\mathrm{dt}\right]+\underset{j=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\λ}}_j{M}_T^{{\mathrm{\α}}_j}$

至此，可以建立求解寿命分布概率密度函数参数的无约束优化模型：

$\left\{\begin{array}{c}\mathrm{Find}:\mathrm{\λ}={[{\mathrm{\λ}}_1,{\mathrm{\λ}}_2,...,{\mathrm{\λ}}_m]}^T\mathrm{and\α}={[{\mathrm{\α}}_1,{\mathrm{\α}}_2,...,{\mathrm{\α}}_m]}^T\\ \mathrm{Minimize}:I(\mathrm{\λ},\mathrm{\α})=\log \left[{\∫}_T\exp (-\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\λ}}_i{t}^{{\mathrm{\α}}_i})\mathrm{dt}\right]+\underset{i=1}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\λ}}_i{M}_T^{{\mathrm{\α}}_i}\end{array}\right.$

上述无约束优化模型的解可以应用拟牛顿法、单纯形法等无约束非线性优化问题求解。得到最优拉格朗日乘子λ＝[λ₀,λ₁,…,λ_m]^T和最优分数矩指数α＝[α₀,α₁,…,α_m]^T。应用MATLAB工具箱中的fminsearch函数进行求解，它的优势在于应用单纯形法而不需要计算目标函数的梯度信息，编程和数值求解方便。

步骤5.5、利用历史故障数据求解核电站反应堆保护系统寿命分布概率密度函数的拉格朗日乘子λ与分数矩指数α；

步骤5.6、将拉格朗日乘子λ＝[λ₀,λ₁,…,λ_m]^T与分数矩指数α＝[α₀,α₁,…,α_m]^T代入步骤5-3中的核电站反应堆保护系统的概率密度函数的解析式估计，得到核电站反应堆保护系统寿命分布概率密度函数，如图5所示。

${\hat{f}}_T\left(t\right)=\exp (-\underset{j=0}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\λ}}_j{t}^{{\mathrm{\α}}_j}),(j=\mathrm{1,2},...,m).$

其中， ${\mathrm{\α}}_0=0,{\mathrm{\λ}}_0=\log \left[{\∫}_T\exp (-\underset{i=0}{\overset{m}{\mathrm{\Σ}}}{\mathrm{\λ}}_i{t}^{{\mathrm{\α}}_i})\mathrm{dt}\right].$

步骤6、利用核电站反应堆保护系统寿命分布概率密度函数求解核电站反应堆保护系统的失效概率，即核电站反应堆保护系统在t时刻之前发生故障的概率(核电站反应堆保护系统寿命分布概率密度函数在时间0到t时刻对t的积分值)，如图6所示。

即： $P_F\left(t\right)=\mathrm{Pr}[T<t]={\hat{F}}_T\left(t\right)={\&Integral;}_0^t{\hat{f}}_X\left(x\right)\mathrm{dx}$

步骤7、利用核电站反应堆保护系统寿命分布概率密度函数求解核电站反应堆保护系统可靠度，即核电站反应堆保护系统在时间t之后仍然正常工作的概率，如图7所示。可靠度与系统的失效概率互补。

即： $R\left(t\right)=\mathrm{Pr}[T\&GreaterEqual;t]=1-{\hat{F}}_T\left(t\right)=1-{\&Integral;}_0^t{\hat{f}}_X\left(x\right)\mathrm{dx}$

步骤8、利用核电站反应堆保护系统寿命分布概率密度函数求解核电站反应堆保护系统动态失效率，即核电站反应堆保护系统工作到时刻t时尚未失效、在时刻t以后的单位时间内发生失效的概率，如图8所示。

即： $\hat{h}\left(t\right)=\frac{{\hat{f}}_T\left(t\right)}{R_T\left(x\right)}=\frac{{\hat{f}}_t\left(t\right)}{1-{\&Integral;}_0^t{\hat{f}}_X\left(x\right)\mathrm{dx}}=\frac{\exp (-\underset{i=0}{\overset{m}{\mathrm{\&Sigma;}}}{\mathrm{\&lambda;}}_i{t}^{{\mathrm{\&alpha;}}_i})}{1-\exp (-\underset{i=0}{\overset{m}{\mathrm{\&Sigma;}}}{\mathrm{\&lambda;}}_i{t}^{{\mathrm{\&alpha;}}_i})},({\mathrm{\&alpha;}}_0=0)$

需要指出的是，图5、图6、图7、图8中MCS表示10⁶次蒙特卡洛数值抽样结果，ME-FM为基于10³个样本的分数矩最大熵优化结果。从中可以看出应用本发明的基于分数矩最大熵优化的可靠性分析方法，仅凭10³个寿命数据样本所得的系统寿命概率分布、可靠性与动态失效率计算精度即与10⁶次蒙特卡洛随机结果精度相同，体现出本发明方法基于小数据样本重构系统寿命分布和动态失效率计算方面的优势与使用价值。

Claims (3)

1.一种核电站反应堆保护系统可靠性分析方法，其特征在于，包括以下步骤：

步骤1、根据核电站反应堆保护系统的结构及功能联系确定其故障树模型；

步骤2、用下行法求解引发顶事件即核电站紧急停堆失效发生的最小割集，即底事件同时发生则导致顶事件发生的组合；

步骤3、获取核电站反应堆保护系统的历史故障数据即各个底事件发生故障的时间统计数据样本；

步骤4、通过最小割集计算核电站反应堆保护系统的寿命统计量，即核电站反应堆保护系统发生故障时已经运行的时间；

步骤5、建立核电站反应堆保护系统寿命分布概率密度函数的信息熵模型，求解核电站反应堆保护系统最优的寿命分布概率密度函数；

步骤6、利用核电站反应堆保护系统寿命分布概率密度函数求解核电站反应堆保护系统的失效概率，即核电站反应堆保护系统在t时刻之前发生故障的概率；

步骤7、利用核电站反应堆保护系统寿命分布概率密度函数求解核电站反应堆保护系统可靠度，即核电站反应堆保护系统在时间t之后仍然正常工作的概率；

步骤8、利用核电站反应堆保护系统寿命分布概率密度函数求解核电站反应堆保护系统动态失效率，即核电站反应堆保护系统工作到时刻t时尚未失效、在时刻t以后的单位时间内发生失效的概率。

2.根据权利要求1所述的核电站反应堆保护系统可靠性分析方法，其特征在于，所述步骤1包括以下步骤：

步骤1.1、确定故障树的顶事件及导致顶事件发生的二级事件；

故障树的顶事件为核电站紧急停堆失效；

导致顶事件的二级事件包括稳压器压力低信号失效、停堆断路器拒开、至少三束控制棒卡住；

任一二级事件发生均导致顶事件发生；

步骤1.2、确定故障树中导致二级事件发生的三级事件；

导致稳压器压力低信号失效事件发生的事件包括稳压器压力传感器失效、三个稳压器压力传感器阈值继电器定值错误；

导致停堆断路器拒开事件发生的事件为两个停堆断路器共因失效；

至少三束控制棒卡住事件视为底事件即导致核电站紧急停堆失效的不能再分的事件；

步骤1.3、确定导致三级事件发生的可能事件即四级事件；

导致稳压器压力传感器失效的事件包括相关的稳压器压力传感器失效；其中任意两个事件同时发生则会导致稳压器压力传感器失效事件的发生；

三个稳压器压力传感器阈值继电器定值错误、两个停堆断路器共因失效均为底事件；

步骤1.4、确定导致四级事件发生的可能事件即五级事件；

步骤1.5、确定导致五级事件发生的可能事件即六级事件，直至该事件不能再分为止，引发顶事件即核电站紧急停堆失效事件发生的所有可能事件搜寻完毕，得到底事件。

3.根据权利要求1所述的核电站反应堆保护系统可靠性分析方法，其特征在于，所述核电站反应堆保护系统寿命分布概率密度函数按如下步骤建立：

步骤5.1、引入最大熵估计方法，建立核电站反应堆保护系统的信息熵模型；

步骤5.2、确定核电站反应堆保护系统信息熵优化模型的约束条件，包括核电站反应堆保护系统寿命统计量的分数矩的建立、核电站反应堆保护系统寿命分布概率密度函数估计的积分值为1；

步骤5.3、引入拉格朗日方程求解核电站反应堆保护系统最大熵约束优化问题，使拉格朗日方程对核电站反应堆保护系统寿命分布概率密度估计量求偏导数，令其值等于0，得到核电站反应堆保护系统寿命分布概率密度函数的解析式估计；

步骤5.4、引入K-L距离方法，建立求解核电站反应堆保护系统寿命分布概率密度函数参数的无约束优化模型；

步骤5.5、利用历史故障数据求解核电站反应堆保护系统寿命分布概率密度函数的拉格朗日乘子λ与分数矩指数α；

步骤5.6、将拉格朗日乘子λ与分数矩指数α代入步骤5-3中的核电站反应堆保护系统的概率密度函数的解析式估计，得到核电站反应堆保护系统寿命分布概率密度函数。