CN112199850A - 一种基于过程特征的系统弹性分布规律研究方法 - Google Patents

Abstract

本发明公开了一种基于过程特征的系统弹性分布规律研究方法，属于系统弹性技术领域。首先确定待研究的系统弹性过程函数及关键参数的概率分布，确定扰动下的弹性度量方法。然后通过预实验确定蒙特卡罗法的仿真次数N₁，并用蒙特卡洛抽样法获取系统弹性过程关键参数的样本，使得样本数量达到仿真次数N₁，再代入待研究的弹性过程函数中，构建出N₁个系统弹性过程，使用牛顿莱布尼茨公式求解弹性积分，根据确定型弹性度量方法计算得到N₁个弹性样本。最后采用直方图对N₁个弹性样本进行理论分布拟合，对拟合结果进行拟合优度检验，直到卡方检验通过。本发明从统计角度重新认识系统弹性属性，根据本发明确定的系统弹性分布规律应用广泛。

Description

一种基于过程特征的系统弹性分布规律研究方法

技术领域

本发明属于系统弹性技术领域，具体是一种基于过程特征的系统弹性分布规律研究方法。

背景技术

弹性(resilience)描述了系统承受扰动以及扰动后的恢复能力，如文献[1]:LiR,Tian X,Yu L,Kang R.A systematic disturbance analysis method for resilienceevaluation:A case study in material handling systems.Sustainability.2019,11(5):1447。在实际应用中，系统不仅会发生各种内部故障，而且还会面临自然灾害或人为攻击等外部干扰。系统对扰动的预测、抵抗、吸收、适应和恢复能力，决定了其性能变化的过程特征。

系统从受到扰动开始至结束的一个完整弹性过程，如图1所示，横轴为时间，纵轴为归一化后的系统性能(为了使不同系统间性能便于比较，将充分运行时的系统性能记为100％，完全失效时记为0％，其他情况按性能百分比记)。图中，Q₀(t)表示系统在正常情况下(未受扰动)的归一化性能，Q(t)表示含扰动行为的归一化系统性能曲线，T_a为用户允许的系统最大恢复时间。t₀时刻系统受到扰动，导致性能Q(t)开始下降，经过时长T_d至t₁时刻，系统性能下降至最低点1-L并同时开始恢复；再经过T_r时长至t₂时刻，系统性能恢复结束，达到新的稳定状态。

系统弹性过程函数以及相关参数分布直接影响了系统弹性表现，如文献[2]：Cimellaro G P,Reinhorn A M,Bruneau M.Framework for analytical quantificationof disaster resilience[J].Engineering Structures.2010,32(11):3639-3649.使用线性、三角和指数三种函数来描述系统的性能下降和恢复过程。亦有众多学者在研究基础设施系统、供应链系统、通风系统或航空系统等系统时，分别总结了其性能降级程度、降级时间和恢复时间等弹性过程关键参数的概率分布规律，如文献[3]：Porter K A,KiremidjianA S,Legrue J S.Assembly-based vulnerability of buildings and its use inperformance evaluation[J].Earthquake Spectra.2001,17(2):291-312.在地震中建筑脆弱性的研究中认为部件的损坏程度服从离散分布；如文献[4]:Li R,Dong Q,Jin C,KangR.a new resilience measure for supply chain networks[J].Sustainability.2017.9(1):1-19.认为供应链系统的节点中断持续时间遵循指数分布，其修复时间通常遵循对数正态分布；文献[5]:Porter K A,Kiremidjian A S,Legrue J S.Assembly-basedvulnerability of buildings and its use in performance evaluation[J].Earthquake Spectra.2001,17(2):291-312.和文献[6]:Myrefelt S.The reliabilityand availability of heating,ventilation and air conditioning systems[J].Energy and Buildings.2004,36(10):1035-1048.在研究工程调度系统和供热、通风与空调系统时发现当故障率为常数时，可修复系统的平均维修时间符合对数正态分布。

由于扰动的随机性，以及系统对扰动响应和恢复过程的随机性，系统弹性是一个随机变量。现有的弹性量化研究主要集中于如何度量弹性，而对于系统弹性这一随机变量的分布规律，却鲜有研究。

发明内容

本发明针对上述问题，给出了一套通过数值仿真分析方法，对已知弹性过程特征函数及参数分布的系统进行弹性分布规律研究的方法，具体是一种基于过程特征的系统弹性分布规律研究方法。

具体步骤如下：

步骤一、针对连续可修复系统，确定待研究的典型弹性过程函数及关键参数的概率分布；

待研究的典型弹性过程函数包括如下三种：

1、线性函数Q_A(t)：性能降级/恢复速度恒定。

式中，t₀表示系统受到扰动并同时开始性能下降的时刻；T_d表示系统性能开始下降到降至最低点的时长；t₁表示系统性能下降至最低点并同时开始恢复的时刻；t₂表示系统性能恢复至扰动发生前状态的时刻；T_r表示系统性能开始恢复到恢复至扰动发生前状态的时长；L表示系统性能下降的程度。

2、三角函数Q_B(t)：性能降级/恢复速度在初始阶段和末尾阶段较慢，中间阶段较快。

3、指数函数Q_C(t)：性能降级速度在前期较快后期较慢，恢复速度反之。

式中，b是指数函数的尺度参数，决定了指数函数偏离线性函数的程度。

弹性分析过程中的关键参数包括：系统性能下降时间T_d、恢复时间T_r以及性能下降程度L。典型概率分布主要有：离散分布、指数分布、威布尔分布、广义极值分布、正态分布和对数正态分布；

步骤二：确定扰动下的系统弹性度量方法；

本发明使用的确定型弹性计算公式如下：

式中，t₀为系统受到扰动并同时开始性能下降的时刻；T_a为用户允许的系统最大恢复时间。Q₀(t)为系统在正常情况下(未受扰动)的归一化性能目标值；Q(t)为弹性过程函数，表示系统在受扰动的情况下随时间变化的归一化性能值；

步骤三：通过预实验确定蒙特卡罗法的仿真次数N₁；

首先，估计仿真次数N₀，进行预实验得到弹性样本

以及弹性样本值的方差S₀ ²；

式中，

为第i次仿真得到的弹性样本值；

为弹性样本的期望，即

然后，利用样本方差S₀ ²，计算满足置信度和仿真误差条件下的最少仿真次数N₁：

式中，取样本方差S₀ ²为方差σ²的无偏估计；

表示标准正态分布的上

分位点，ε为仿真误差，1-α为置信度；

步骤四：使用蒙特卡洛抽样法获取系统弹性过程关键参数的样本，使得样本数量达到仿真次数N₁；

根据关键参数的概率分布进行蒙特卡洛抽样，得到N₁组抽样样本：T_d(i)，T_r(i)，L(i)(i＝1,2,...,N₁)。

其中，样本的常规分布通过已有的随机数生成器生成，包括：指数分布随机数生成器exprnd(λ)、正态分布随机数生成器normrnd(μ,σ)或均匀分布随机数生成器unifrnd(a,b)等。

对于未提供随机数生成器的分布类型，通过反函数法进行抽样，具体为：首先通过随机数生成器生成一组符合0-1均匀分布的随机数Y_K，再根据待抽样随机变量的累计分布函数F(t)，通过反函数方式生成符合其分布的随机数样本，即X_K＝F^-1(Y_K)。

步骤五、将N₁组关键参数样本，代入待研究的弹性过程函数中，构建出N₁个系统弹性过程；并使用牛顿莱布尼茨公式求解弹性积分，根据确定型弹性度量方法计算N₁个弹性样本

步骤六：采用直方图对N₁个弹性样本

进行理论分布拟合；

首先，根据拇指规则来确定样本分组数目，如下：

k＝[1+3.3lgN₁] (7)

式中，k为分组数目；

然后，构建确定型弹性度量值

的样本分布直方图，通过观察其特征，初步识别出几种理论分布：如果其形状呈对称或近似对称分布，那么样本可能符合正态分布或形状参数在3-4之间的威布尔分布；如果数据是右偏的，那么可能符合指数分布、对数正态分布或威布尔分布等；如果数据是分散的，那么可能符合均匀分布或离散分布。

最后，通过概率图进一步验证，选择坐标点最接近直线的分布类型，使用极大似然估计法进行参数估计。

步骤七：使用卡方检验方法对拟合结果进行拟合优度检验，判断是否通过检验，如果是，则得出结论；否则，重复步骤六直到卡方检验通过。

本发明的优点与积极效果在于：

(1)本发明一种基于过程特征的系统弹性分布规律研究方法，考虑了扰动的随机性以及系统受扰动后性能降级和恢复行为的随机性，根据弹性过程特征，设计了系统弹性分布规律研究方法，有助于从统计角度重新认识系统弹性属性；

(2)本发明一种基于过程特征的系统弹性分布规律研究方法，确定出的系统弹性分布规律可用于：①为弹性参数选取提供参考，②为系统弹性分析提供新思路，③为弹性测试方案制定提供依据。

附图说明

图1是现有技术中采用的系统弹性过程示意图；

图2是本发明一种基于过程特征的系统弹性分布规律研究方法的流程图；

图3是本发明确定的待研究的典型弹性过程函数示意图；

图4是本发明实施例中弹性样本的经验分布直方图；

图5是本发明实施例中弹性样本拟合广义极值分布和威布尔分布的概率密度图；

图6是本发明实施例中弹性样本拟合广义极值分布和威布尔分布的概率图。

具体实施方式

下面将结合附图和实施例对本发明作进一步的详细说明。

本发明提供了一种基于过程特征的系统弹性分布规律研究方法，通过蒙特卡洛抽样获取系统弹性过程关键参数，进而计算系统弹性，最后通过理论分布拟合来确定这种情况下的系统弹性分布规律。如图2所示，具体步骤如下：

步骤一、针对连续可修复系统，确定待研究的典型弹性过程函数及关键参数的概率分布；

在研究基于过程特征的系统弹性分布规律前，首先要确定系统受到扰动后的性能降级和恢复过程函数，以及描述弹性过程的关键参数(即系统性能下降时间T_d、恢复时间T_r以及性能下降程度L)的概率分布特征。

下面分别给出典型的弹性过程函数和弹性过程关键参数概率分布以作参考。

本发明根据文献[2]提出的弹性过程函数，结合本发明的参数表达，给出三种典型弹性过程函数，如图3所示，如下：

1、线性函数Q_A(t)：性能降级/恢复速度恒定。作为性能降级/恢复函数，适用于缺乏系统响应、应急准备、可用资源等相关信息的系统。

式中，t₀表示系统受到扰动并同时开始性能下降的时刻；T_d表示系统性能开始下降到降至最低点的时长；t₁表示系统性能下降至最低点并同时开始恢复的时刻；t₂表示系统性能恢复至扰动发生前状态的时刻；T_r表示系统性能开始恢复到恢复至扰动发生前状态的时长；L表示系统性能下降的程度。

2、三角函数Q_B(t)：性能降级/恢复速度在初始阶段和末尾阶段较慢，中间阶段较快。作为性能降级函数，使用该函数的系统通常是扰动在初始阶段对系统影响较小，但随着作用时间/强度增加，扰动影响变大，末期又由于系统设计部分性能不易受扰动影响，使得其性能降级速度变缓；作为性能恢复函数，使用该函数的系统通常是恢复初始阶段资源缺乏或有限，随着资源配置到位优先恢复对系统性能影响大的部分，后恢复对系统性能影响小的部分，因此刚开始恢复速度能迅速提升，末期又会略微放缓。

3、指数函数Q_C(t)：性能降级速度在前期较快后期较慢，恢复速度反之。作为性能降级函数，适用于扰动开始时对系统性能影响较小，但随着扰动作用时间/强度的增加，对系统影响迅速增大的系统；作为性能恢复函数，适用于初始资源迅速注入，优先恢复对系统影响大的部分导致前期恢复速度较快，后恢复影响小的部分导致后期速度较缓的系统。

式中，b是指数函数的尺度参数，决定了指数函数偏离线性函数的程度。

弹性分析过程中的关键参数包括：系统性能下降时间T_d、恢复时间T_r以及性能下降程度L。典型概率分布根据不同系统的特征确定，主要有以下几种：

(1)离散分布，该类型样本的分布律如下：

P(x＝i)＝p_i,i＝1,2,...,n (11)

式中，p_i为样本取值为i的概率；

(2)指数分布，该类型样本的经验分布函数如下：

式中，λ为特征参数；

(3)威布尔分布，该类型样本的经验分布函数如下：

式中，β为形状参数，θ为尺度参数；

(4)广义极值分布，该类型样本的分布函数如下：

式中，ξ为形状参数，σ为尺度参数，μ为位置参数；

(5)正态分布，该类型样本的概率密度函数如下：

式中，μ为数学期望，σ为方差；

(6)对数正态分布，该类型样本的概率密度函数如下：

式中，μ为数学期望，σ为方差；

步骤二：确定扰动下的系统弹性度量方法；

本发明使用的确定型弹性计算公式如下：

式中，Q₀(t)为系统在正常情况下(未受扰动)的归一化性能目标值；由于没有扰动时系统水平可能是在百分之百附近波动的，所以这里也用函数来表示。Q(t)为弹性过程函数，表示系统在受扰动的情况下的归一化性能值；t₀为系统受到扰动并开始性能下降的时刻；T_a为用户允许的系统最大恢复时间，该参数反映了系统受扰动后在要求恢复时间内的平均性能水平满足目标要求的程度，适用于单次扰动下的弹性计算。

步骤三：通过预实验确定蒙特卡罗法的仿真次数N₁；

首先，估计一个比较合理的仿真次数N₀(如1000次)，进行预实验得到弹性样本：

以及弹性样本值的方差S₀ ²；

式中，

为第i次仿真得到的弹性样本值；

为弹性样本均值，即

然后，假定要求在置信度水平为1-α时，仿真误差为±ε，即

由于使用蒙特卡洛法从一个总体中抽样获得的样本X是独立且同分布的随机变量。根据中心极限定理，可推导出满足置信度和仿真误差条件下的最少仿真次数N₁，计算公式如下：

式中，取样本方差S₀ ²为方差σ²的无偏估计；

表示标准正态分布的上

分位点；

步骤四、使用蒙特卡洛抽样获取系统弹性过程中包括性能下降时间T_d、恢复时间T_r以及性能下降程度L的样本，使得样本数量达到仿真次数N₁；

根据三个关键参数的概率分布进行蒙特卡洛抽样，得到一组抽样样本：T_d(i)，T_r(i)，L(i)(i＝1,2,...,N₁)。

其中，样本的常规分布可以通过已有的随机数生成器生成，如MATLAB中提供了指数分布随机数生成器exprnd(λ)、正态分布随机数生成器normrnd(μ,σ)、均匀分布随机数生成器unifrnd(a,b)等。

对于未提供随机数生成器的分布类型，通过反函数法进行抽样，具体为：首先通过随机数生成器生成一组符合0-1均匀分布的随机数Y_K，再根据待抽样随机变量的累计分布函数F(t)，通过反函数方式生成符合其分布的随机数样本，即X_K＝F^-1(Y_K)；

步骤五、将关键参数的N₁组样本，代入待研究的弹性过程函数中，构建出N₁个系统弹性过程；并使用牛顿莱布尼茨公式求解弹性积分，根据确定型弹性度量方法计算N₁个弹性样本

步骤六：采用直方图对N₁个弹性样本

进行理论分布拟合；

首先，根据拇指规则来确定样本分组数目，如下：

k＝[1+3.3lgN₁] (20)

式中，k为分组数目；[]代表取整，N₁为样本总量。

然后，构建确定型弹性度量值

的样本分布直方图，通过观察其特征，初步识别可能的理论分布：如果其形状呈对称或近似对称分布，那么样本可能符合正态分布或形状参数在3-4之间的威布尔分布；如果数据是右偏的，那么可能符合指数分布、对数正态分布或威布尔分布等；如果数据是分散的，那么可能符合均匀分布或离散分布。

最后，通过概率图进一步验证，选择拟合程度最好的分布类型(即坐标点最接近直线的一种分布)，使用极大似然估计法进行参数估计。

步骤七：使用卡方检验方法对拟合结果进行拟合优度检验，判断是否通过检验，如果是，则得出结论；否则，重复步骤六直到卡方检验通过。

卡方检验步骤如下：

(1)作出两种假设：

原假设H₀：弹性样本

服从A分布；备择假设H₁：弹性样本

不服从A分布。

(2)根据样本分组数k计算组距：通过样本中最大值与最小值的差除以分组数k统计每组的频数。

(3)根据A分布的分布函数，利用F_A(x_j)-F_A(x_j-1)求出样本落入区间[x_j-1,x_j]的概率p_j。

(4)为了表示样本观测值与理论分布之间的偏离程度，使用下式计算得到卡方检验值χ²：

式中，O_j为第j组中的样本量即频数；p_j为如果假设H₀成立，样本

落入第j组的概率。

(5)根据给定的显著性水平α，由卡方分布相应分位点求出临界值

其中m为估计参数的个数。然后，判断该临界值是否大于χ²，如果是，则接受原假设H₀，得出结论：弹性样本

服从A分布；否则接受备择假设H₁，重新进行理论分布拟合直至通过卡方检验。

实施例：

本发明以下面这类系统为例说明本发明的应用过程：

首先，系统弹性过程函数为线性函数；三个关键参数的选择如下：系统性能下降时间T_d和恢复时间T_r分别遵循特征参数为60的指数分布；系统性能下降程度L遵循离散均匀分布，分布律为：P(L＝0.01i)＝0.01,i＝1,2,…,100。

已知Q₀(t)＝1，用户允许的系统最大恢复时间T_a＝800s，并要求在置信度1-α＝95％的前提下仿真误差不超过ε＝0.1％。

本实施例估计一个比较合理的仿真次数N₀＝1000，进行预实验得到1000个弹性样本，然后计算样本方差S₀ ²＝0.1296，查表得到

根据仿真次数的公式计算得到N₁≥4979，取仿真次数N₁＝5000。

接着借助MATLAB中的随机数生成函数exprnd(60)和0.01unidrnd(100)分别对性能下降时间T_d、恢复时间T_r以及性能下降程度L进行5000次随机抽样，获取5000组随机变量T_d,T_r,L。然后分别将这5000组数据代入Q_A(t)获得扰动后的系统性能函数，并利用确定型弹性度量值的公式(17)计算得到5000个系统弹性样本。

接下来，利用公式k＝[1+3.3lgN₁]计算出样本分组数k＝13，观察其样本分布直方图，如图4所示，发现其可能符合广义极值分布或威布尔分布。

借助MATLAB的分布拟合工具箱进行拟合。其中，广义极值分布拟合得到的形状参数ξ为-1.0261，尺度参数σ为0.0420，位置参数μ为0.9591；威布尔分布拟合得到的尺度参数θ为0.9774，形状参数β为41.2611。分析还得到样本拟合广义极值分布和威布尔分布的概率密度对比图，如图5所示；以及如图6所示的概率图。从两图中皆可看出拟合广义极值分布的曲线与样本点重合程度较高，初步说明广义极值分布拟合效果最好。

假设H₀：弹性样本

遵循广义极值分布；H₁：弹性样本

遵循广义极值分布。

合并样本量小于5的分组，统计每组的频数O_j。根据广义极值分布的分布函数F_gev(x)，计算F_gev(x_j)-F_gev(x_j-1)求出样本落入区间[x_j-1,x_j]的概率p_j，计算得而出卡方检验值。卡方检验过程如表1所示。

表1

本实施例中自由度为k-m-1＝10-3-1＝6，查表可得显著性水平α＝0.05时卡方临界值为12.6,χ²＝11.3589小于该临界值，于是接受假设H₀：系统弹性

符合广义极值分布。

Claims (4)

1.一种基于过程特征的系统弹性分布规律研究方法，其特征在于，具体步骤如下：

步骤一、针对连续可修复系统，确定待研究的典型弹性过程函数及关键参数的概率分布；

待研究的典型弹性过程函数包括如下三种：

1、线性函数Q_A(t)：性能降级/恢复速度恒定；

式中，t₀表示系统受到扰动并同时开始性能下降的时刻；T_d表示系统性能开始下降到降至最低点的时长；t₁表示系统性能下降至最低点并同时开始恢复的时刻；t₂表示系统性能恢复至扰动发生前状态的时刻；T_r表示系统性能开始恢复到恢复至扰动发生前状态的时长；L表示系统性能下降的程度；

2、三角函数Q_B(t)：性能降级/恢复速度在初始阶段和末尾阶段较慢，中间阶段较快；

3、指数函数Q_C(t)：性能降级速度在前期较快后期较慢，恢复速度反之；

式中，b是指数函数的尺度参数，决定了指数函数偏离线性函数的程度；

弹性分析过程中的关键参数包括：系统性能下降时间T_d、恢复时间T_r以及性能下降程度L；典型概率分布主要有：离散分布、指数分布、威布尔分布、广义极值分布、正态分布和对数正态分布；

步骤二：确定扰动下的系统弹性度量方法；

本发明使用的确定型弹性计算公式如下：

式中，t₀为系统受到扰动并同时开始性能下降的时刻；T_a为用户允许的系统最大恢复时间；Q₀(t)为系统在正常情况下(未受扰动)的归一化性能目标值；Q(t)为弹性过程函数，表示系统在受扰动的情况下随时间变化的归一化性能值；

步骤三：通过预实验确定蒙特卡罗法的仿真次数N₁；

步骤四：使用蒙特卡洛抽样法获取系统弹性过程关键参数的样本，使得样本数量达到仿真次数N₁；

步骤五、将N₁组关键参数样本，代入待研究的弹性过程函数中，构建出N₁个系统弹性过程；并使用牛顿莱布尼茨公式求解弹性积分，根据确定型弹性度量方法计算N₁个弹性样本

步骤六：采用直方图对N₁个弹性样本

进行理论分布拟合；

步骤七：使用卡方检验方法对拟合结果进行拟合优度检验，判断是否通过检验，如果是，则得出结论；否则，重复步骤六直到卡方检验通过。

2.如权利要求1所述的一种基于过程特征的系统弹性分布规律研究方法，其特征在于，所述的步骤三具体为：

首先，估计仿真次数N₀，进行预实验得到弹性样本

以及弹性样本值的方差S₀ ²；

式中，

为第i次仿真得到的弹性样本值；

为弹性样本的期望，即

然后，利用样本方差S₀ ²，计算满足置信度和仿真误差条件下的最少仿真次数N₁：

式中，取样本方差S₀ ²为方差σ²的无偏估计；

表示标准正态分布的上

分位点，ε为仿真误差，1-α为置信度。

3.如权利要求1所述的一种基于过程特征的系统弹性分布规律研究方法，其特征在于，所述的步骤四为：

根据关键参数的概率分布进行蒙特卡洛抽样，得到N₁组抽样样本：T_d(i)，T_r(i)，L(i)(i＝1,2,...,N₁)；

其中，样本的常规分布通过已有的随机数生成器生成，包括：指数分布随机数生成器exprnd(λ)、正态分布随机数生成器normrnd(μ,σ)或均匀分布随机数生成器unifrnd(a,b)；

对于未提供随机数生成器的分布类型，通过反函数法进行抽样，具体为：首先通过随机数生成器生成一组符合0-1均匀分布的随机数Y_K，再根据待抽样随机变量的累计分布函数F(t)，通过反函数方式生成符合其分布的随机数样本，即X_K＝F^-1(Y_K)。

4.如权利要求1所述的一种基于过程特征的系统弹性分布规律研究方法，其特征在于，所述的步骤六具体如下：

首先，根据拇指规则来来确定样本分组数目，如下：

k＝[1+3.3lgN₁] (7)

式中，k为分组数目；

然后，构建确定型弹性度量值

的样本分布直方图；通过观察其特征，初步识别出几种理论分布：如果其形状呈对称或近似对称分布，那么样本可能符合正态分布或形状参数在3-4之间的威布尔分布；如果数据是右偏的，那么可能符合指数分布、对数正态分布或威布尔分布等；如果数据是分散的，那么可能符合均匀分布或离散分布；

最后，通过概率图进一步验证，选择坐标点最接近直线的分布类型，使用极大似然估计法进行参数估计。

Images