CN107729612A - 一种印刷装备关键部件的隐患评估方法 - Google Patents

Abstract

本发明提供一种印刷装备关键部件的隐患评估方法，首先通过计算滚动轴承振动加速度数据的状态特征变量的正常域边界，确定其运行状态正常域和异常域空间，然后基于正常域边界和安全裕度进行滚动轴承的安全评估。本发明能够准确地识别印刷装备关键部件是否发生异常，为印刷装备的状态监测、风险评价、故障诊断提供了有效的方法。

Description

一种印刷装备关键部件的隐患评估方法

技术领域

本发明涉及印刷装备状态监测技术领域，尤其涉及一种印刷装备关键部件的隐患评估方法。

背景技术

目前，多数的大型印刷装备尚未应用在线状态识别和监控技术，仍是采用基于人工的离线和定期检查手段，仅仅能够做到事故后处理和定性分析，无法对装备和其中关键零部件进行实时、动态、定量的运行状态评估。此外，由于印后装备具有多样化和更新快的特点，传统的检查手段无法跟上装备功能和结构的升级速度。因此，建立一套相对普适的大型印刷装备隐患状态评估和分析方法是保障印刷生产的安全高效和最大化发挥装备效能的必要工作。

要保障大型印刷装备整体系统的安全运行，必须首要保证其中关键和核心零部件处于正常运行状态。据统计，滚动轴承、齿轮箱等关键部件是大型印刷装备在内的复杂机械装备中使用最频繁且故障率最高的零部件，因此，从关键部件的状态监测和隐患评估入手，首先完成运行安全关键零部件的隐患评估，是当前亟待解决的首要关键技术问题。

目前，国内外针对印刷装备的隐患评估研究较少，相关研究工作更多地集中在故障处理和信息采集等方面。在印刷装备关键部件的状态评估方面，在实际的生产现场，工程师们主要根据印刷品质量问题(如墨杠、磨毛等)通过专家经验对印刷装备旋转部件的工作状态进行粗判；在理论研究方面，主要从印品质量的角度关注旋转部件对印刷压力和纸张张力等的影响，现有成果基本集中在基于动力学分析和基于振动测试的直接方法，以及基于图像识别的间接方法。

目前，针对“隐患”的研究成果主要集中两个方面：一是在矿业和电气等领域的事故隐患分级、风险评估体系、企业安全管理策略研究等宏观方面；二是针对某类设备运行状态的早期故障辨识和隐患监测等微观问题。目前，后者相关研究成果较为丰富，而前者在印刷装备领域未见开展。

发明内容

本发明的目的是针对印刷装备关键部件隐患评估问题，提出了基于区域估计的方法。本发明可以为印刷装备关键部件的隐患辨识和评估提供一种新方法，为解决印刷装备和系统的实时状态评估提供技术支持和指导。

本发明将具体研究对象(如滚动轴承或齿轮箱等某一印刷装备关键部件)的运行状态分为正常、异常两种情况。用运行状态变量所构成的一个空间表征研究对象的全部运行状态，而其中的正常和异常两种情况分别对应整个空间的不同区域，正常状态对应的空间区域认定为正常域，异常状态对应的空间区域认定为异常域。

本发明定义正常域(Normal Region,NR)是一种从域的角度描述系统整体可正常稳定运行区域的定量模型。在印刷装备关键部件状态评估的研究中，正常域是一个在研究对象各运行状态相关变量(如振动、速度、位移、提取的状态特征等)所确定的空间内，用于评价研究对象运行状态是否正常的区域。具体到某一具体的研究对象，正常域是指在该对象运行状态相关变量(不同的研究对象，其运行状态相关变量亦不相同)所确定的空间内，用于评价当前对象的运行状态是否正常的区域。

若对象的运行状态落在正常域内，则表示对象的运行状态是健康的正常的，隐患等级较低，若当对象的运行状态落在除正常域以外的其他区域，可能有两种情况，其一是对象已经处于故障的运行状态，且这一故障将导致非安全事件(如事故)的发生，此时隐患等级高，其二是对象虽然已经处于不健康不正常的运行状态，但其不会直接导致非安全事件的发生，此时隐患等级中等。

如图1所示，若研究对象运行状态位于正常域(图中绿色所示区域)内，则可判定研究对象运行是正常的，否则认为研究对象运行异常；若有隐患发生或外部干扰，可能导致对象运行状态临近正常域边界，说明对象可能存在隐患，有发生故障的趋势，需要及时采取相关预防控制措施，使对象运行状态远离正常域边界(如图中虚线箭头所示)；若对象运行状态在正常域边界的临近位置时，没有及时正确地采取防控措施，则对象运行状态将会进一步恶化，最终发生故障，演化至异常域(如图中故障轨迹所示)，进而导致事故发生。

从区域与空间划分的角度来说，是通过估计空间中的边界，由边界将空间划分成不同区域。当仅需区分正常状态和异常状态时，只需估计出一条正常域边界(如图1)，将运行状态相关变量空间划分为两个区域即可。此时的区域边界可用如下的分类决策函数表示即可

f(X)＝sign[Bound(X)] (1)

其中：X＝(x₁,x₂,…,x_n)∈Rⁿ，为运行状态相关变量空间中的状态特征变量，n为状态特征变量的维数，x₁,x₂,…,x_n为各个状态特征变量的元素值；Bound(X)为用于划分两个不同区域的边界函数，即用于划分两个区域的边界方程可用Bound(X)＝0来描述。

本发明中的安全裕度，将其定义为：对于某一具体研究对象，在其运行状态相关变量空间内，由某时刻该对象各运行状态相关变量取值所确定的状态点距离正常域边界的最小距离。其中，关于距离这一概念，需要指出的是，为避免某一确定的距离概念或计算方法可能不适用于某些特定印刷装备的在线状态评价，并为给出较普适的安全裕度的概念，故在此定义中的“距离”指广义距离，对于不同的研究对象可选择合适的距离计算方法，如欧式距离、马氏距离等。此外，正常域边界在二维运行状态相关变量空间内为一曲线，在三维运行状态相关变量空间内为一曲面，在更高维的运行状态相关变量空间内为可将其理解为一超平面，则安全裕度则可理解为状态点距离曲线/曲面/超平面的最小距离。

本发明是在明确上述正常域、异常域、正常域边界、安全裕度等相关概念的基础上，提出隐患评估步骤。针对某一具体的研究对象，从正常域估计的角度出发，完成其隐患评估大致分为两步。

第一步：计算研究对象的正常域边界，确定运行状态正常域和异常域空间

正常域边界由正常域边界方程Bound(X)＝Bound(x₁,x₂,…,x_n)＝c来描述，其中x₁,x₂,…,x_n为可表征对象动态行为的多个内部变量，n为变量个数，Bound(X)为表征对象运行状态的输出变量(假设其与正常度量成正相关，即正常程度越高Bound(X)的值越大)，c为表征正常阈值的一常数。当n＝3时，正常域边界方程Bound(x₁,x₂,x₃)＝c即所描述的正常域边界为三维空间中的某一曲面，如图2所示。当n＝1时，此正常域边界方程即退化为简单的单变量阈值。定义Bound(X)>c所表示的区域(假设正常域边界方程所确定曲面以下的区域)为正常域；反之Bound(X)<c所表示的区域(假设正常域边界方程所确定曲面以上的区域)为异常域。

第二步：基于正常域边界和安全裕度进行对象的安全评估

基于已经确定的正常域边界、正常域和异常域空间，首先由对象运行状态变量的数据计算Bound(X)的值，依据Bound(X)是否小于c判断对象当前运行状态是否位于正常域内；其次，若对象运行状态处于正常域内，则计算对象运行状态点距离正常域边界的距离，此距离的大小定义为安全裕度，安全裕度越大表明对象运行状态的正常程度越高，若对象运行状态处于异常域内，则立即给出报警信息；最后，给出定量化的评估结果。

具体采用如下技术方案：

(1)获取滚动轴承的振动加速度数据作为原始状态数据；

(2)将所述原始状态数据进行分段后再进行经验模式分解，每一段数据分解后对应一个分解矩阵，将分解矩阵进行多元统计分析，计算每个分解矩阵的T²统计量控制限和SPE统计量控制限两个数值，将此两个数值作为每段振动加速度数据的状态特征变量；

(3)估计正常域边界：将所有状态特征变量进行正常域的边界估计，正常域边界由正常域边界方程Bound(X)＝Bound(x₁,x₂,…,x_n)＝c来描述，其中x₁,x₂,…,x_n为可表征对象动态行为的多个内部变量，n为变量个数，Bound(X)为表征对象运行状态的输出变量，c为表征正常阈值的常数；

(4)由滚动轴承状态变量计算Bound(X)的值，依据Bound(X)是否小于c判断对象当前运行状态是否位于正常域内；若滚动轴承运行状态处于正常域内，将滚动轴承的状态特征变量点到正常域边界的欧式距离作为安全裕度，在某时刻的运行状态相关变量空间内的运行状态点为A(a₁,a₂,…,a_m)，同一空间中的正常域边界上任一点为B(b₁,b₂,…,b_m)，最小欧式距离问题的最优化模型为：

采用离散牛顿法的方程组求解上述模型,计算出当前状态特征变量点距离正常域边界的欧式距离，得到安全裕度，安全裕度越大表明滚动轴承运行状态的正常程度越高，若滚动轴承运行状态处于异常域内，则立即给出报警信息。

优选地，离散牛顿法的方程组求解模型的具体方式为：

对于给定点Q，若在隐式曲线上P点达到局部极大或极小距离时，则向量QP与P点的切向量T_P垂直，对于给定点Q，若在隐式曲面上某点M处达到局部极大距离或局部极小距离时，则有QM向量与M点处的法向N_M是平行的，将二维平面隐式曲线表示为

f(x,y)＝0

T_P为隐式曲线上点P的切向量，则

T_P＝(-f_y,f_x)

若给定点Q到隐式曲线上的P点达到局部极值距离，则有

Q_P*T_P＝0

找出点Q到隐式曲线上能达到极值距离的点P₁,P₂,…,P_k，求解给定点Q到隐式曲线的最小距离d_min，其中点P₁,P₂,…,P_k满足非线性方程组：

牛顿法求解上述非线性方程组的计算过程如下：

步骤1：设隐式曲线C的端点分别为E₁和E₂，取初始值X＝[x₁,x₂,…,x_n]，t，h>0，0<t<1，计算f_i(X)→E₁(i)，其中i＝1,2,…,n；

步骤2：若有则方程组的一组实数解为X＝[x₁,x₂,…,x_n]^T，计算过程结束，否则继续步骤3；

步骤3：计算f_i(X_j)→E₂(i,j)，i,j＝1,2,…,n，其中X_j＝[x₁,…x_j-1,x_j+h,x_j+1,…,x_n]；

步骤4：解线性代数方程组E₂Z＝E₁，其中Z＝[z₁,z₂,…,z_n]，并计算

步骤5：计算x_i-hz_i/β→x_i，i＝1,2,…,n；

步骤6：t×h→h，转步骤1。

本发明具有如下有益的技术效果：

本发明提供了一种基于区域估计的印刷装备关键部件隐患评估方法，能够准确地识别印刷装备关键部件是否发生异常，为印刷装备的状态监测、风险评价、故障诊断提供了有效的方法。

附图说明

下面结合附图和实例对本发明进一步说明：

图1是区域估计示意图。

图2是基于正常域边界估计的隐患评估示意图。

图3是隐患评估过程。

图4是安全裕度计算示例。

具体实施方式

下面结合附图对本发明的实施例具体说明，本发明中以滚动轴承这一印刷装备中故障率发生率最高的关键部件为实例进行展开，其中常用的轴承型号包括SKF 51113、SKF33216等。如图3所示，印刷装备关键部件隐患评估方法包括以下具体实施步骤：

第一：针对具体的研究对象，掌握其基本原理和工况，分析其特点，如是否便于建模、是否能够获取状态监测数据等情况。研究状态特征的表达，充分考虑现场的工作环境和对象特点，选取便于获取同时能够充分地灵敏地反映对象运行状态变化的特征量。本实施例中针对滚动轴承，考虑其基本工况，获取其振动加速度数据作为原始状态数据。

第二：依据所选取的状态特征变量，获取对象的运行状态特征，计算得到状态特征变量值，即得到运行状态相关变量空间中的状态点，必要时需要考虑噪声及不确定环境对研究对象及获取过程的影响。本实施例中将滚动轴承的振动数据进行分段后再进行经验模式分解，每一段数据分解后对应一个分解矩阵，将分解矩阵进行多元统计分析，计算每个分解矩阵的T²统计量控制限和SPE统计量控制限两个数值，此两个数值作为每段振动加速度数据的特征变量。

第三：针对实际的运行状态数据积累情况，完成仅有正常状态情况下的单值正常域边界估计或具有完备的正常和异常多状态数据的多值正常域边界估计，将运行状态相关变量空间划分为所需的正常和异常两个区域。本实施例中单值正常域边界估计采用典型的如凸包生成、支持向量数据描述等单值分类方法，多值正常域边界估计采用典型的支持向量机或C均值聚类等智能分类方法，所得的分类决策函数即为正常域边界估计。

第四：根据已完成的边界估计，将状态点对应至某一区域，计算状态点至正常域边界的广义距离作为安全裕度。本实施例中，将滚动轴承的状态特征变量点到正常域边界的欧式距离作为安全裕度，采用基于几何特征的方程求解方法。设某装备在某时刻的运行状态相关变量空间内的运行状态点为A(a₁,a₂,…,a_m)，同一空间中的正常域边界上任一点为B(b₁,b₂,…,b_m)，则最小欧式距离问题的求解可建立如下最优化模型：

如式(1)中的约束条件所示，可将正常域边界方程看作一个隐式曲线/曲面C，以下基于计算复杂度较低的离散牛顿法的方程组求解方法，以点到隐式曲线的最小距离求解为例，给出常用的具体算法。

对于给定点Q，若在隐式曲线上P点达到局部极大或极小距离时，则向量QP与P点的切向量T_P垂直，可得一个方程组。对于给定点Q，若在隐式曲面上某点M处达到局部极大距离或局部极小距离时，则应有QM向量与M点处的法向N_M是平行的，因此得另一个方程组。将二维平面曲线隐式表示为

f(x,y)＝0 (2)

T_P为隐式曲线(2)上点P的切向量，则

T_P＝(-f_y,f_x) (3)

若给定点Q到隐式曲线(2)上的P点达到局部极值距离，则有

Q_P*T_P＝0 (4)

可见，求解给定点Q到隐式曲线(2)的最小距离d_min的关键在于找到点Q到曲线(2)上能达到极值距离的点P₁,P₂,…,P_k。因这些点同时满足曲线方程(2)和方程(4)，故可得非线性方程组：

下面讨论方程组(5)的解法，给出常用的离散牛顿法求解方法。设有非线性方程组

其中f_i(x₁,x₂,...,x_n)，i＝1,2,…,n，为实变量非线性函数，是已确定的多元函数。上式可用向量的形式表示为

在求解式(7)时常用的迭代方法是牛顿法，其处理过程如下。

对于非线性方程组F(X)＝0，其中F(X)＝[f₁(x),f₂(x),…,f_n(x)]^T，f_i(x)的偏导数矩阵记为J(X)或F’(X)。

设X^*为F(X)＝0的解，设X^(k)＝[x₁ ^(k),x₂ ^(k),…,x_n ^(k)]^T为X^*的近似解，则利用多元函数f_i(x)在X^(k)的泰勒公式，有

其中C_i在X^(k)于X所连的线段内。

若采用式(9)中的线性函数P_i(X)近似代替f_i(x)，并将线性方程组

的解作为X^*的第k+1此近似解X^(k+1)。将式(10)转换为矩阵形式。

F(X^(k))+J(X^(k))(X-X^(k))＝0 (11)

若J(X^(k))是非奇异矩阵，则可得牛顿迭代公式

采用牛顿法求解非线性方程组(7)，可用如下形式。

故，由式(13)可知，每次计算X^(k)→X^(k+1)时，需要三步：计算矩阵J(X^(k))和F(X^(k))；求解线性方程组J(X^(k))ΔX^(k)＝F(X^(k))；计算X^(k+1)＝X^(k)-ΔX^(k)。但在实际问题中，很多情况下J(X^(k))的计算相当复杂，因此可用相应的差商代替J(X^(k))的元素，即：

其中，h足够小，且

则J(X^(k))ΔX^(k)＝F(X^(k))变换为

经式(14)至式(17)的变换，牛顿法求解非线性方程组的计算过程如下。

步骤1：设隐式曲线C的端点分别为E₁和E₂，取初始值X＝[x₁,x₂,…,x_n]，t，h>0，0<t<1，计算f_i(X)→E₁(i)，其中i＝1,2,…,n；

步骤2：若有则方程组的一组实数解为X＝[x₁,x₂,…,x_n]^T，计算过程结束，否则继续步骤3；

步骤3：计算f_i(X_j)→E₂(i,j)，i,j＝1,2,…,n，其中X_j＝[x₁,…x_j-1,x_j+h,x_j+1,…,x_n]；

步骤4：解线性代数方程组E₂Z＝E₁，其中Z＝[z₁,z₂,…,z_n]，并计算

步骤5：计算x_i-hz_i/β→x_i，i＝1,2,…,n；

步骤6：t×h→h，转步骤1。

以下是本发明的一个实施例。

实施例中所需数据是由Dr.Kenneth A.Loparo提供的滚动轴承实验数据，轴承型号为SKF205-2RS型深沟球轴承，电机负载3马力，转速1730r/min(约28.8r/s)，振动加速度数据采集点事负载端，采样频率48k Hz，共采样4次，每次采样时间为10s；采样时间10s。

本实施例采用采样数据按照前述具体实施步骤，对原始数据进行分段和经验模式分解，提取每段数据分解后的T²统计量控制限和SPE统计量控制限两个值作为状态特征，获得2维平面空间上的状态点，如图4中实心点所示。其中的正常域边界已经基于此轴承的正常和异常状态数据采用支持向量机分类的方法获得，如图4所示中黑色实曲线所示。计算当前状态点距离正常域边界的最小距离，即可获得安全裕度。

经计算当前运行状态的T²统计量控制限和SPE统计量控制限分别为X₁＝26.86和X₂＝147.01，即当前状态点在正常域平面的坐标为(26.86,147.01)。然后使用前述的基于几何特征的方程求解方法，计算状态点(26.86,147.01)与正常域边界曲线的最小欧式距离，找到该状态点距离曲线上点(21.40,139.40)的欧式距离最小，为9.358，则安全裕度的值为9.358。

进一步地，计算另个时刻的状态点与正常域边界的相对位置，经计算第二个状态点的T²统计量控制限和SPE统计量控制限分别为X₁＝25.60和X₂＝87.27，即其在正常域平面的坐标为(25.60，87.27)。仍使用基于几何特征的方程求解方法计算该点与正常域边界曲线的最小距离。可找到点(25.60，87.27)距离曲线上点(24.10，89.10)的欧式距离最小，为2.366，即此时安全裕度的值为2.366。可见，同一对象在不同时刻的安全裕度分别为9.358和2.366，说明后一时刻的状态与前一时刻的状态相比，其安全性程度变差，状态点有向异常域转移的趋势，隐患风险加大。

Claims (2)

1.一种印刷装备关键部件的隐患评估方法，其特征在于，包括如下步骤：

(1)获取滚动轴承的振动加速度数据作为原始状态数据；

(2)将所述原始状态数据进行分段后再进行经验模式分解，每一段数据分解后对应一个分解矩阵，将分解矩阵进行多元统计分析，计算每个分解矩阵的T²统计量控制限和SPE统计量控制限两个数值，将此两个数值作为每段振动加速度数据的状态特征变量；

(3)估计正常域边界：将所有状态特征变量进行正常域的边界估计，正常域边界由正常域边界方程Bound(X)＝Bound(x₁,x₂,…,x_n)＝c来描述，其中x₁,x₂,…,x_n为可表征对象动态行为的多个内部变量，n为变量个数，Bound(X)为表征对象运行状态的输出变量，c为表征正常阈值的常数；

(4)由滚动轴承状态变量计算Bound(X)的值，依据Bound(X)是否小于c判断对象当前运行状态是否位于正常域内；若滚动轴承运行状态处于正常域内，将滚动轴承的状态特征变量点到正常域边界的欧式距离作为安全裕度，在某时刻的运行状态相关变量空间内的运行状态点为A(a₁,a₂,…,a_m)，同一空间中的正常域边界上任一点为B(b₁,b₂,…,b_m)，最小欧式距离问题的最优化模型为：

<mfenced open = "{" close = ""> <mtable> <mtr> <mtd> <mi>min</mi> </mtd> <mtd> <mrow> <mi>d</mi> <mi>i</mi> <mi>s</mi> <mi>tan</mi> <mi>c</mi> <mi>e</mi> <mrow> <mo>(</mo> <mi>A</mi> <mo>,</mo> <mi>B</mi> <mo>)</mo> </mrow> <mo>=</mo> <msqrt> <mrow> <munderover> <mo>&amp;Sigma;</mo> <mrow> <mi>i</mi> <mo>=</mo> <mn>1</mn> </mrow> <mi>m</mi> </munderover> <msup> <mrow> <mo>(</mo> <msub> <mi>a</mi> <mi>i</mi> </msub> <mo>-</mo> <msub> <mi>b</mi> <mi>i</mi> </msub> <mo>)</mo> </mrow> <mn>2</mn> </msup> </mrow> </msqrt> </mrow> </mtd> </mtr> <mtr> <mtd> <mrow> <mi>s</mi> <mi>t</mi> <mo>.</mo> </mrow> </mtd> <mtd> <mrow> <mi>B</mi> <mi>o</mi> <mi>u</mi> <mi>n</mi> <mi>d</mi> <mrow> <mo>(</mo> <mi>B</mi> <mo>)</mo> </mrow> <mo>-</mo> <mi>c</mi> <mo>=</mo> <mi>B</mi> <mi>o</mi> <mi>u</mi> <mi>n</mi> <mi>d</mi> <mrow> <mo>(</mo> <msub> <mi>b</mi> <mn>1</mn> </msub> <mo>,</mo> <msub> <mi>b</mi> <mn>2</mn> </msub> <mo>,</mo> <mn>...</mn> <mo>,</mo> <msub> <mi>b</mi> <mi>m</mi> </msub> <mo>)</mo> </mrow> <mo>-</mo> <mi>c</mi> <mo>=</mo> <mn>0</mn> </mrow> </mtd> </mtr> </mtable> </mfenced>

采用离散牛顿法的方程组求解上述模型,计算出当前状态特征变量点距离正常域边界的欧式距离，得到安全裕度，安全裕度越大表明滚动轴承运行状态的正常程度越高，若滚动轴承运行状态处于异常域内，则立即给出报警信息。

2.如权利要求2所述的印刷装备关键部件的隐患评估方法，其特征在于，离散牛顿法的方程组求解模型的具体方式为：

对于给定点Q，若在隐式曲线上P点达到局部极大或极小距离时，则向量QP与P点的切向量T_P垂直，对于给定点Q，若在隐式曲面上某点M处达到局部极大距离或局部极小距离时，则有QM向量与M点处的法向N_M是平行的，将二维平面隐式曲线表示为

f(x,y)＝0

T_P为隐式曲线上点P的切向量，则

T_P＝(-f_y,f_x)

若给定点Q到隐式曲线上的P点达到局部极值距离，则有

Q_P*T_P＝0

找出点Q到隐式曲线上能达到极值距离的点P₁,P₂,…,P_k，求解给定点Q到隐式曲线的最小距离d_min，其中点P₁,P₂,…,P_k满足非线性方程组：

<mfenced open = "{" close = ""> <mtable> <mtr> <mtd> <mi>f</mi> <mo>(</mo> <mi>x</mi> <mo>,</mo> <mi>y</mi> <mo>)</mo> <mo>=</mo> <mn>0</mn> </mtd> </mtr> <mtr> <mtd> <msub> <mi>Q</mi> <mi>P</mi> </msub> <mo>*</mo> <msub> <mi>T</mi> <mi>P</mi> </msub> <mo>=</mo> <mn>0</mn> </mtd> </mtr> </mtable> </mfenced>

牛顿法求解上述非线性方程组的计算过程如下：

步骤1：设隐式曲线C的端点分别为E₁和E₂，取初始值X＝[x₁,x₂,…,x_n]，t，h>0，0<t<1，计算f_i(X)→E₁(i)，其中i＝1,2,…,n；

步骤2：若有则方程组的一组实数解为X＝[x₁,x₂,…,x_n]^T，计算过程结束，否则继续步骤3；

步骤3：计算f_i(X_j)→E₂(i,j)，i,j＝1,2,…,n，其中X_j＝[x₁,…x_j-1,x_j+h,x_j+1,…,x_n]；

步骤4：解线性代数方程组E₂Z＝E₁，其中Z＝[z₁,z₂,…,z_n]，并计算

步骤5：计算x_i-hz_i/β→x_i，i＝1,2,…,n；

步骤6：t×h→h，转步骤1。