CN113552804B - 一种基于qr-frols算法的转子系统narx建模方法 - Google Patents

Abstract

本发明提供一种基于QR‑FROLS算法的转子系统NARX建模方法，包括：采集用于表征转子系统动力学特性的振动响应信号，并将其作为建立NARX模型的输入和输出；基于输入和输出的时滞变量，构建模型项候选集；对模型项候选集中的每一个模型项对应列做归一化处理，得到模型项矩阵；对归一化处理后的模型项矩阵做列主元QR分解；基于R矩阵对角线值的累计占比，对模型项进行筛选，形成新的模型项候选集；基于新的模型项候选集，采用FROLS算法建立NARX模型。本发明的技术方案能够有效降低转子系统谐波信号时滞变量组成的模型项候选集的多重共线性，进而解决了利用FROLS算法对转子系统谐波信号建模结果不稳定以及预测输出易发散的问题。

Description

一种基于QR-FROLS算法的转子系统NARX建模方法

技术领域

本发明涉及模型构建技术领域，具体而言，尤其涉及一种基于QR-FROLS算法的转子系统NARX建模方法。

背景技术

转子系统的动力学分析离不开合理的动力学模型，通过机理建模方式建立复杂转子系统模型时，需要考虑转动惯量、陀螺力矩、载荷等多种模型参数。模型精度取决于是否对转子动力学影响因素考虑充足，以及对实际转子系统中的基本规律的全面把握。数据驱动建模方法只根据研究对象中的输入输出数据之间的相互关系直接建模，无需专业知识以及复杂的计算过程，模型能适用于高度非线性和严重不确定性系统，从而为解决转子系统建模问题提供了一条有效途径。

NARX(Nonlinear Auto Regressive with eXogenous inputs非线性有源自回归)模型建立了当前输出与过去输入、输出之间的关系，其具体的定义如下：

y(k)＝F[y(k-1),y(k-2),...,y(k-n_y),u(k-1),u(k-2),...,u(k-n_u)]+e(k)

其中，y(k)、u(k)和e(k)分别代表系统的输出输入和噪声序列；n_y和n_u分别代表系统的最大延迟；F[·]表示一个非线性函数，针对不同的系统，一般可以选择不同类型的非线性函数，现在普遍使用的都是多项式函数。利用系统输入、输出数据建立NARX模型常用的算法是FROLS(Forward Regression Orthogonal Least Squares向前回归正交最小二乘)算法，其基本原理是通过ERR(Error Reduction Ratio误差减少比)指标从模型项候选集中逐次选取出对输出影响最重要的模型项。NARX模型因为其结构简单，普适性好，鲁棒性强，构建容易，能反映真实系统特性等特点，在诸如冶金，机械，医疗等的许多领域得到了广泛应用。

现有的利用NARX模型表征研究对象动力学特性的研究大部分是以随机信号作为系统激励，因为随机信号包含不同的幅值和频率特性。如专利《基于非线性输出频率响应函数的结构损伤检测方法》(CN102520070A)中利用经功率放大器放大后的白噪声信号驱动激振器激励桥梁模型结构，并利用加速度传感器测量桥梁结构在白噪声激励下的输出信号，以此建立桥梁结构NARX模型。但由于随机激励并不适用于转子系统，文献《The NARXModel-Based System Identification on Nonlinear,Rotor-Bearing Systems》(DOI:10.3390/app7090911)中研究者通过采集转子系统升速过程中产生的多谐波信号建立NARX模型表征转子系统来解决上述问题，所得NARX模型可以表征转子系统谐波信号频率范围内的时域及频域特征。

但是，针对上述现有技术，还存在如下缺陷：

转子系统的振动信号多为谐波信号，跟随机信号相比，谐波信号时滞变量组成的模型项之间具有较高的相关性，使得回归建模过程中自变量对因变量的影响部分重叠或存在包含作用，这种现象称为多重共线性。逐步回归方法可以在一定程度上筛选存在多重共线性的自变量组合中对因变量解释较大的变量，但这种方法在自变量共线性较为严重时不能完全解决问题，甚至对自变量个数以及模型长度十分敏感。所以当转子系统谐波信号的模型项之间的相关性较强时，使用FROLS算法可能会造成辨识结果中出现冗余模型项，进而导致模型预测输出发散的现象。

发明内容

根据上述提出的技术问题，提供一种基于QR-FROLS算法的转子系统NARX建模方法。本发明首先通过列主元QR分解对候选模型项矩阵进行预处理；再通过处理结果筛选出具有代表性的模型项组成新的候选集；最后利用FROLS算法筛选出对系统输出贡献最大的前若干个模型项组成转子系统NARX模型。经仿真及实验验证，本发明方法能够有效降低转子系统谐波信号时滞变量组成的模型项候选集的多重共线性，进而解决了利用FROLS算法对转子系统谐波信号建模结果不稳定，预测输出易发散的问题。

本发明采用的技术手段如下：

一种基于QR-FROLS算法的转子系统NARX建模方法，包括：

采集用于表征转子系统动力学特性的振动响应信号，并将其作为建立NARX模型的输入和输出；

基于所述输入和输出的时滞变量，构建模型项候选集；

对模型项候选集中的每一个模型项对应列做归一化处理，得到模型项矩阵；

对归一化处理后的模型项矩阵做列主元QR分解；

基于R矩阵对角线值的累计占比，对模型项进行筛选，形成新的模型项候选集；

基于新的模型项候选集，采用FROLS算法建立NARX模型。

进一步地，所述采集用于表征转子系统动力学特性的振动响应信号，并将其作为建立NARX模型的输入和输出，包括：

采集振动响应信号的转速范围为[ω_min,ω_max]；

采集振动响应信号的采样频率f_s为大于被采样信号中所关心最高频率的2.56倍。

进一步地，所述基于所述输入和输出的时滞变量，构建模型项候选集，包括：

选取建立NARX模型的输入时滞阶数n_u和输出时滞阶数n_y，得到输入输出时滞变量：

[y(k-1),y(k-2),...,y(k-n_y),u(k-1),u(k-2),...,u(k-n_u)]

其中，y表示输出，u表示输入，k表示时间；

选取模型非线性阶数l，利用上述时滞变量生成非线性阶数小于等于l的模型项p_m(t)，其中m∈[1,M]，M是生成的模型项总数：

M＝[(n_y+n_u)+l]！/[(n_y+n_u)！l！]

生成模型项候选集{p_m(t),m∈[1,M]}。

进一步地，所述对模型项候选集中的每一个模型项对应列做归一化处理，得到模型项矩阵，包括：

将所有模型项合成模型项矩阵，如下：

对模型项矩阵的每一列进行归一化处理，将所有模型项向量中的元素映射到[-1,1]之间，避免列主元QR分解过程受到模型项向量模值的影响：

上式中，i表示第i个时刻的模型项序列，也是模型项矩阵P的第i行，m表示第m个模型项，也代表模型项矩阵P的第m列；P_(i,m)表示归一化之前模型项矩阵P的第(i,m)个元素，

表示归一化处理后的模型项矩阵的第(i,m)个元素，其上角标n表示英文单词“normalized”。

进一步地，所述对归一化处理后的模型项矩阵做列主元QR分解，包括：

将归一化处理后的模型项矩阵Pⁿ分解成如下形式：

Pⁿp＝QR

其中，Q表示分解后的正交矩阵，R表示分解后的上三角矩阵，p表示列变换索引。

进一步地，所述基于R矩阵对角线值的累计占比，对模型项进行筛选，形成新的模型项候选集，包括：

当模型项矩阵P列主元QR分解后的R矩阵对角线元素diag(R)从大到小排列时，将diag(R)第一次累计占总体比例达到99.9％时的前J列对应的模型项保留作为新的模型项候选集：

其中，上角标r表示英文单词“reserved”。

进一步地，所述基于新的模型项候选集，采用FROLS算法建立NARX模型，包括：

建立NARX模型包括模型的结构选择和参数估计；

模型的结构选择即确定一种模型结构，用于准确描述系统输入输出变量之间的函数关系；

参数估计即在模型结构确定的基础上，计算模型结构对应的系数。

进一步地，所述采用FROLS算法建立NARX模型，包括：

在新的模型项候选集

中，首先引入最重要的模型项，然后通过添加次重要的模型项来细化模型，确定模型的结构，估计模型中的未知参数。

进一步地，所述方法还包括对采用FROLS算法建立的NARX模型进行验证的步骤，包括：

将所述NARX模型用MPO(Model Predicted Output)验证方法进行模型验证，系统的MPO模型验证输出为：

其中，y表示预测输出，

表示实际输出，u表示实际输入，a、b、c表示系数；相同输入下NARX模型验证输出与实际系统输出之间的差异可以量化NARX模型的拟合效果进而判断NARX模型对系统的表征能力。

计算模型验证指标NMSE，其表达式如下：

其中，y_real(t)表示系统真实输出，y(t)表示NARX模型预测输出，N表示数据点数。

较现有技术相比，本发明具有以下优点：

1、本发明提供的基于QR-FROLS算法的转子系统NARX建模方法，通过对模型项矩阵做列主元QR分解，进而对候选模型项进行筛选，同时保证了候选模型项集合的完备性和精简性，解决了转子系统振动谐波信号时滞变量组成的模型项之间相关性较大的问题，使得最终建立的NARX模型中模型项既对系统输出有统计学影响又相互独立，从而解决了预测输出易发散的问题。

2、本发明提供的基于QR-FROLS算法的转子系统NARX建模方法，相比于传统FROLS算法，本发明的方法特点在于降低了时滞阶数和模型项长度对建模结果收敛性的影响，增强了算法的鲁棒性，更适合在先验条件不足，转子系统特征未知的情况下应用。

基于上述理由本发明可在模型构建等领域广泛推广。

附图说明

为了更清楚地说明本发明实施例或现有技术中的技术方案，下面将对实施例或现有技术描述中所需要使用的附图做以简单地介绍，显而易见地，下面描述中的附图是本发明的一些实施例，对于本领域普通技术人员来讲，在不付出创造性劳动性的前提下，还可以根据这些附图获得其他的附图。

图1为本发明方法流程图。

图2为本发明实施例提供的输入信号拼接结果曲线图。

图3为本发明实施例提供的输出信号拼接结果曲线图。

图4为本发明实施例提供的NARX模型拟合效果图。

图5为本发明实施例提供的diag(R)变化规律示意图。

图6为本发明实施例提供的NARX模型拟合效果图。

具体实施方式

为了使本技术领域的人员更好地理解本发明方案，下面将结合本发明实施例中的附图，对本发明实施例中的技术方案进行清楚、完整地描述，显然，所描述的实施例仅仅是本发明一部分的实施例，而不是全部的实施例。基于本发明中的实施例，本领域普通技术人员在没有做出创造性劳动前提下所获得的所有其他实施例，都应当属于本发明保护的范围。

需要说明的是，本发明的说明书和权利要求书及上述附图中的术语“第一”、“第二”等是用于区别类似的对象，而不必用于描述特定的顺序或先后次序。应该理解这样使用的数据在适当情况下可以互换，以便这里描述的本发明的实施例能够以除了在这里图示或描述的那些以外的顺序实施。此外，术语“包括”和“具有”以及他们的任何变形，意图在于覆盖不排他的包含，例如，包含了一系列步骤或单元的过程、方法、系统、产品或设备不必限于清楚地列出的那些步骤或单元，而是可包括没有清楚地列出的或对于这些过程、方法、产品或设备固有的其它步骤或单元。

如图1所示，本发明提供了一种基于QR-FROLS算法的转子系统NARX建模方法，包括：

S1、采集用于表征转子系统动力学特性的振动响应信号，并将其作为建立NARX模型的输入和输出；

S2、基于所述输入和输出的时滞变量，构建模型项候选集；

S3、对模型项候选集中的每一个模型项对应列做归一化处理，得到模型项矩阵；

S4、对归一化处理后的模型项矩阵做列主元QR分解；

S5、基于R矩阵对角线值的累计占比，对模型项进行筛选，形成新的模型项候选集；

S6、基于新的模型项候选集，采用FROLS算法建立NARX模型。

将本发明方法应用到碰磨故障下的转子试验系统，通过盘左和盘右两处的振动响应表征转子系统，研究碰磨故障下的转子系统特征。将两个电涡流传感器布置在转盘一左一右两个位置的水平方向上，连接LMS数据采集系统对转子实验系统进行数据采集。

具体实施时，作为本发明优选的实施方式，所述步骤S1的具体实现过程，包括：

采集振动响应信号的转速范围为[ω_min,ω_max]；

采集振动响应信号的采样频率f_s为大于被采样信号中所关心最高频率的2.56倍。

在本实施例中，在200r/min、300r/min、400r/min和500r/min四个转速下采集碰磨转子系统盘左和盘右两个位置的振动响应信号，采样频率为2048Hz。

对每个位置下测得的响应数据进行拼接，按照转速由小到大的顺序，每个转速下的数据取前五千个数据点拼接成一段数据。以盘左的拼接数据作为输入信号记为u(t)，如图2所示。以盘右的拼接数据作为输出信号记为y(t)，如图3所示。所得NARX模型可以近似表征转子系统频率在[200,500]r/min内的频域特征。

具体实施时，作为本发明优选的实施方式，所述步骤S2的具体实现过程，包括：

选取建立NARX模型的输入时滞阶数n_u和输出时滞阶数n_y，得到输入输出时滞变量：

[y(k-1),y(k-2),...,y(k-n_y),u(k-1),u(k-2),...,u(k-n_u)]

其中，y表示输出，u表示输入，k表示时间；

选取模型非线性阶数l，利用上述时滞变量生成非线性阶数小于等于l的模型项p_m(t)，其中m∈[1,M]，M是生成的模型项总数：

M＝[(n_y+n_u)+l]！/[(n_y+n_u)！l！]

生成模型项候选集{p_m(t),m∈[1,M]}。

在本实施例中，选取建立NARX模型的输入输出时滞阶数n_u＝5和n_y＝5，得到输入输出时滞变量：

选取模型非线性阶数l＝3，利用上述时滞变量生成模型项M＝285个，采用全部模型项作为候选模型项集{p_m(t),m∈[1,285]}。设置预期模型项个数为13个，直接采用FROLS算法建立NARX模型，结果如下：

y(t)＝2.509y(t-1)-0.106u(t-5)+0.214y(t-5)-0.236u(t-1)-10.687y²(t-1)y(t-5)+0.699y(t-1)y²(t-3)+7.980y²(t-2)y(t-5)-43.434y³(t-1)-1.721y(t-2)+76.580y²(t-1)y(t-2)+0.341u(t-2)-33.195y³(t-2)+2.010y(t-1)y²(t-5)

NARX模型的MPO验证输出与系统真实输出的拟合效果如图4所示。其预测输出发散，建模失败。

具体实施时，作为本发明优选的实施方式，所述步骤S3的具体实现过程，包括：

将所有模型项合成模型项矩阵，如下：

对模型项矩阵的每一列进行归一化处理，将所有模型项向量中的元素映射到[-1,1]之间，避免列主元QR分解过程受到模型项向量模值的影响：

上式中，i表示第i个时刻的模型项序列，也是模型项矩阵P的第i行，m表示第m个模型项，也代表模型项矩阵P的第m列；P_(i,m)表示归一化之前模型项矩阵P的第(i,m)个元素，

表示归一化处理后的模型项矩阵的第(i,m)个元素，其上角标n表示英文单词“normalized”。

具体实施时，作为本发明优选的实施方式，所述步骤S4的具体实现过程，包括：

将归一化处理后的模型项矩阵Pⁿ分解成如下形式：

Pⁿp＝QR

其中，Q表示分解后的正交矩阵，R表示分解后的上三角矩阵，p表示列变换索引。三角矩阵R对角线值diag(R)的大小变化规律如图5所示。

QR分解是由Gram–Schmidt正交化推理出来的一种方阵分解形式，列主元QR分解中经过列变换后，三角矩阵R对角线值diag(R)的大小会自动由大到小排列，列变换索引p也代表着原模型项矩阵中对应的模型项序列与其他序列相关性呈由小到大的程度排列。

具体实施时，作为本发明优选的实施方式，所述步骤S5的具体实现过程，包括：

本发明中，QR-FROLS算法第一次从候选模型项集{p_m(t),m∈[1,M]}中筛选模型项的依据是模型项在候选集中是否具有代表性，同时保证下一步的FROLS算法在完备的候选集中搜寻模型项。出于对候选模型项集的完备性和精简性两者之间平衡，本发明中，当模型项矩阵P列主元QR分解后的R矩阵对角线元素diag(R)从大到小排列时，将diag(R)第一次累计占总体比例达到99.9％时的前J列对应的模型项保留作为新的模型项候选集：

其中，上角标r表示英文单词“reserved”。

模型项相关性较大的情况下，一小部分模型项就能代表全部模型项后选集，累计占比会十分快速达到99.9％。模型项相关性较小的情况，累计占比达到99.9％时会保留大部分模型项，但因为这种条件下建模不易出现发散的情况，所以即使截断较少的模型项依然能够保持良好的建模效果。所以此方法的有效性不受模型项保留个数J的影响，只与R矩阵的对角线值diag(R)的累计占比大小有关。

本实施例中，将diag(R)第一次累计占总体比例达到99.9％时的前J＝47个模型项保留作为筛选后的模型项候选集

含有的模型项如下：

具体实施时，作为本发明优选的实施方式，所述步骤S6的具体实现过程，包括：

建立NARX模型包括模型的结构选择和参数估计；在新的模型项候选集

中，首先引入最重要的模型项，然后通过添加次重要的模型项来细化模型，确定模型的结构，估计模型中的未知参数。

模型的结构选择即确定一种模型结构，用于准确描述系统输入输出变量之间的函数关系；

参数估计即在模型结构确定的基础上，计算模型结构对应的系数。

具体实施时，作为本发明优选的实施方式，所述步骤S6之后还包括步骤S7，对采用FROLS算法建立的NARX模型进行验证，包括：

模型验证是系统辨识研究的基础部分之一，用以检验数学模型是否为实际系统的无偏估计。将所述NARX模型用MPO(Model Predicted Output)验证方法进行模型验证，系统的MPO模型验证输出为：

其中，y表示预测输出，

表示实际输出，u表示实际输入，a、b、c表示系数；相同输入下NARX模型验证输出与实际系统输出之间的差异可以量化NARX模型的拟合效果进而判断NARX模型对系统的表征能力。在MPO模型验证中，下一步的输出序列，是由上一步的模型验证输出得到，每步的误差会逐渐累积，是一种较为严格的模型验证方法。

计算模型验证指标NMSE(Normalized Mean Square Error，归一化均方误差)，NMSE是一种较为常用的模型误差量化指标，其表达式如下：

其中，y_real(t)表示系统真实输出，y(t)表示NARX模型预测输出，N表示数据点数。

采用实际工况下的数据建模结果中，NMSE指标小于0.1一般认为模型有能力表征对象系统。仿真信号由于没有噪音以及各种干扰成分的存在，建模结果中NMSE指标要求应比实际工况数据更严格。

本实施例中，以

为候选集，同样设置预期模型项个数为13个，采用FROLS算法建立转子系统NARX模型，结果如下：

y(t)＝1.158y(t-1)-0.041u(t-5)-0.145y(t-5)-0.013u(t-1)-16.294y²(t-2)y(t-5)+13.281y³(t-3)-0.170u(t-3)+0.219u(t-2)+15.639u(t-1)y(t-1)y(t-5)-1.074u²(t-1)y(t-3)-12.508u(t-2)y²(t-2)-3.051u(t-1)y²(t-5)+3.374u(t-5)y²(t-1)

NARX模型的MPO验证输出与系统真实输出的拟合效果如图6所示。

模型验证指标NMSE计算结果如下所示，说明本发明方法对转子系统NARX建模预测结果收敛性有明显改善效果。

最后应说明的是：以上各实施例仅用以说明本发明的技术方案，而非对其限制；尽管参照前述各实施例对本发明进行了详细的说明，本领域的普通技术人员应当理解：其依然可以对前述各实施例所记载的技术方案进行修改，或者对其中部分或者全部技术特征进行等同替换；而这些修改或者替换，并不使相应技术方案的本质脱离本发明各实施例技术方案的范围。

Claims (6)

1.一种基于QR-FROLS算法的转子系统NARX建模方法，其特征在于，包括：

采集用于表征转子系统动力学特性的振动响应信号，并将其作为建立NARX模型的输入和输出；

基于所述输入和输出的时滞变量，构建模型项候选集，包括：

选取建立NARX模型的输入时滞阶数n_u和输出时滞阶数n_y，得到输入输出时滞变量：

[y(k-1),y(k-2),...,y(k-n_y),u(k-1),u(k-2),...,u(k-n_u)]

其中，y表示输出，u表示输入，k表示时间；

选取模型非线性阶数l，利用上述时滞变量生成非线性阶数小于等于l的模型项p_m(t)，其中m∈[1,M]，M是生成的模型项总数：

M＝[(n_y+n_u)+l]！/[(n_y+n_u)！l！]

生成模型项候选集{p_m(t),m∈[1,M]}；

对模型项候选集中的每一个模型项对应列做归一化处理，得到模型项矩阵，包括：

将所有模型项合成模型项矩阵，如下：

对模型项矩阵的每一列进行归一化处理，将所有模型项向量中的元素映射到[-1,1]之间，避免列主元QR分解过程受到模型项向量模值的影响：

上式中，i表示第i个时刻的模型项序列，也是模型项矩阵P的第i行，m表示第m个模型项，也代表模型项矩阵P的第m列；P_(i,m)表示归一化之前模型项矩阵P的第(i,m)个元素，

表示归一化处理后的模型项矩阵的第(i,m)个元素，其上角标n表示英文单词“normalized”；

对归一化处理后的模型项矩阵做列主元QR分解，包括：

将归一化处理后的模型项矩阵Pⁿ分解成如下形式：

Pⁿp＝QR

其中，Q表示分解后的正交矩阵，R表示分解后的上三角矩阵，p表示列变换索引；

基于R矩阵对角线值的累计占比，对模型项进行筛选，形成新的模型项候选集；

基于新的模型项候选集，采用FROLS算法建立NARX模型。

2.根据权利要求1所述的基于QR-FROLS算法的转子系统NARX建模方法，其特征在于，所述采集用于表征转子系统动力学特性的振动响应信号，并将其作为建立NARX模型的输入和输出，包括：

采集振动响应信号的转速范围为[ω_min,ω_max]；

采集振动响应信号的采样频率f_s为大于被采样信号中所关心最高频率的2.56倍。

3.根据权利要求1所述的基于QR-FROLS算法的转子系统NARX建模方法，其特征在于，所述基于R矩阵对角线值的累计占比，对模型项进行筛选，形成新的模型项候选集，包括：

当模型项矩阵P列主元QR分解后的R矩阵对角线元素diag(R)从大到小排列时，将diag(R)第一次累计占总体比例达到99.9％时的前J列对应的模型项保留作为新的模型项候选集：

其中，上角标r表示英文单词“reserved”。

4.根据权利要求1所述的基于QR-FROLS算法的转子系统NARX建模方法，其特征在于，所述基于新的模型项候选集，采用FROLS算法建立NARX模型，包括：

建立NARX模型包括模型的结构选择和参数估计；

模型的结构选择即确定一种模型结构，用于准确描述系统输入输出变量之间的函数关系；

参数估计即在模型结构确定的基础上，计算模型结构对应的系数。

5.根据权利要求4所述的基于QR-FROLS算法的转子系统NARX建模方法，其特征在于，所述采用FROLS算法建立NARX模型，包括：

在新的模型项候选集

中，首先引入最重要的模型项，然后通过添加次重要的模型项来细化模型，确定模型的结构，估计模型中的未知参数。

6.根据权利要求1所述的基于QR-FROLS算法的转子系统NARX建模方法，其特征在于，所述方法还包括对采用FROLS算法建立的NARX模型进行验证的步骤，包括：

将所述NARX模型用MPO(Model Predicted Output)验证方法进行模型验证，系统的MPO模型验证输出为：

其中，y表示预测输出，

表示实际输出，u表示实际输入，a、b、c表示系数；相同输入下NARX模型验证输出与实际系统输出之间的差异可以量化NARX模型的拟合效果进而判断NARX模型对系统的表征能力；

计算模型验证指标NMSE，其表达式如下：

其中，y_real(t)表示系统真实输出，y(t)表示NARX模型预测输出，N表示数据点数。

Images