CN110826184B - 一种在时变时滞下narx模型结构和参数的变分贝叶斯辨识方法 - Google Patents

Abstract

本发明公开了一种针对时变时滞下NARX模型的结构和参数的辨识方法。该辨识方法主要包括：首先采用加权多项式作为基函数来表达NARX模型，利用稀疏估计的基础，通过向参数先验中引入稀疏因子独立作用于每个子模型的权重，在整个变分贝叶斯框架下正确地将模型的结构选择出来。对于模型参数的辨识，将每一时刻的时滞值作为缺失变量，利用变分贝叶斯的迭代公式，随着模型结构的正确选择，将对应的未知参数和时变时滞估计出来。本发明的优点：(1)在时变时滞下，能将NARX模型的结构和参数有效地辨识出来(2)能实现对每一时刻时滞的分布进行估计。

Description

一种在时变时滞下NARX模型结构和参数的变分贝叶斯辨识 方法

技术领域

本发明属于系统辨识领域，涉及带有时变时滞的非线性自回归模型的辨识方案

背景技术

近年来，非线性系统辨识成为一个重要而具有挑战性的问题。具有外部输入的非线性自回归(NARX)模型是一类建立在带外部输入的线性自回归模型的基础上的非线性黑箱模型。NARX可以描述一般的非线性，具有良好的函数逼近能力，因此其辨识受到广泛关注。NARX可通过子模型加权的方式表示为参数线性的紧凑形式，子模型可选择为多项式、径向基函数或小波函数。因此NARX模型的辨识包括选择一个具有好的解释能力的最简模型结构及其参数辨识。

目前，NARX模型的识别方法可分为几种类型。其中，最小二乘法和最大似然法是常用的经典方法，包括正回归法、正、反修剪法、稀疏估计法和EM算法。此外，贝叶斯模型识别在近几十年发展迅速，它比其他算法有许多优势。例如：1.不确定性被内部描述，对于分析、仿真和控制设计都很有用2.过度拟合可通过对过于复杂的模型进行自然惩罚来避免3.即使对于样本相对较少的数据记录，模型不确定性也可以精确地量化4.可以在可用的地方加入信息。 NARX识别的贝叶斯方法包括基于高斯过程的非参数方法和基于可逆跳马尔可夫链蒙特卡罗(RJMCMC)等。近年来，一些学者将稀疏技术与贝叶斯方法相结合，发展了稀疏贝叶斯学习方法、稀疏贝叶斯增广拉格朗日算法等。此外， W.R.Jacobs等人还提出了一种新的稀疏贝叶斯识别方法，该方法采用稀疏诱导超先验和变分推理，比MCMC方法快一个数量级。

现提出的NARX模型的辨识方法大都假设没有时滞存在或者只存在单一时滞，但是由于网络传输、硬件限制、化学反应过程等原因，许多控制过程存在时变时滞。如何在时变时滞下对NARX模型结构和参数辨识，成为亟待解决的技术问题。对此，D.H.Zhou等人提出了一种改进的强跟踪滤波器(MSTF)来估计 NARX过程的时变延迟和参数。然而延时变化率必须比输入信号慢，否则算法将失效。期望最大化(EM)算法是识别具有时滞等潜在变量的系统的有效方法。变分贝叶斯(VB)作为EM的一种推广，可以实现参数和时间延迟的估计分布，而不是点估计。因此，在时变延迟下，有学者提出利用VB算法对ARX线性系统进行辨识，可以取得较好的效果,同时也对将线性系统推广到非线性系统提供了很好的思路。

发明内容

对于在时变时滞下NARX模型的结构和参数辨识问题，为了克服现有技术的不足，本发明提供了一种辨识方案。本发明的目的是通过以下技术方案实现的：对带输入时变时滞的NARX模型，使用带稀疏因子的变分贝叶斯方法来辨识结构和参数。首先采用加权多项式作为基函数来表达NARX模型，利用ARD稀疏估计的基础，通过向参数先验中引入稀疏因子独立作用于每个子模型的权重，在整个变分贝叶斯框架下正确地将模型的结构选择出来。对于模型参数的辨识，将每一时刻的时滞值作为隐变量，利用变分贝叶斯的迭代公式，随着模型结构的正确选择，将对应的未知参数和时变时滞估计出来。

本发明的流程图如图1所示，其特征在于包括两大阶段：

第一阶段，进行结构预设/结构初始化，将NARX非线性模型转化为多项式子模型加权和的形式；

第一阶段所述的带时变时滞的单输入单输出NARX模型描述如下：

其中，f(·)是某个非线性函数；{u_k}和{y_k}分别是采样得到的输入和输出数据；下标k表示第k个采样时刻；n_u,n_y分别是输入和输出的最大动态阶次；λ_k为第k个采样时刻出现的时滞，

{v_k}是均值为0，方差为δ^-1的高斯分布白噪声；

可对动态阶次n_u、n_y及多项式阶次n_l设置上限，将NARX模型表示为多项式子模型加权和的形式：

Φ_k＝[φ₁,φ₂,…,φ_M]，θ＝[θ₁,θ₂,…,θ_M]^T

其中，φ_m是由

构成的多项式子模型，例如，

M表示多项式子模型的个数由输入输出的动态阶次n_u、n_y以及多项式的最大阶次n_l决定，θ_m表示每个子模型的权重参数；因此可得一个初始的模型集Φ_k。初始的最大模型集中会出现冗余子模型结构项，需将其从子模型集中剔除。

第二阶段，对多项式子模型加权和的形式进行结构辨识和参数辨识，即选出能够表达系统的最佳多项式子模型集Φ_k，以及求出时变时滞下模型对应的参数 {Θ,β}，其中Θ＝{θ,δ,α}，β＝{β₁,β₂,…,β_j,…}，θ表示子模型的权重向量，δ表示噪声方差的逆，α表示参数θ的先验分布中引入的稀疏因子，β_j表示时滞值j 出现的概率，具体包括以下步骤：

步骤一：设置参数的先验分布，所述的设置参数的先验分布指对使用变分贝叶斯方法进行辨识时所需要的参数设置参数先验，包括设置参数θ的先验分布，α的先验分布，噪声方差的逆δ的先验分布，以及时滞的先验分布，具体如下：

为实现结构辨识，在参数θ的先验分布中引入稀疏因子α，同时为实现共轭分布，选择θ的先验分布为正态分布，具体如下，

p(θ|α)＝p(θ|0,(A)^-1)

其中，α＝(α₁,α₂,…,α_m,…,α_M)^T＝diag(A)，α_m用于控制θ_m变化的幅度，由分析可知，如果

则θ_m越接近于0，说明第m个子模型项对数据来说不显著，为冗余结构，应该将其从模型集中移除；

同样，为得到一个共轭分布，将α的先验分布设置为Gamma分布：

考虑到辨识问题的一般情况，噪声的方差未知，同样将方差的逆δ作为需要辨识的参数并设置先验分布，即如下：

p(δ)＝Gamma(c₀,d₀) 假设时滞出现的最大值为D，设置时滞的先验分布，可将每种时滞出现的概率初始化为相同的值β_j即如下：

p(λ_k＝j)＝β_j

其中，{a₀,b₀,c₀,d₀,D}是需要人为给定的初始化超参数；

步骤二：在VB框架下，在步骤一的基础上进行特定结构下的参数辨识，再用其中的稀疏因子对多项式子模型加权和的形式进行结构修剪即把冗余的多项式子模型剔除，用修剪后的新结构再次进行参数辨识，直到多项式子模型加权和的形式最终只保留一项子模型；计算每种结构和及其对应参数下的Evidence Lower Bound，即F[q(λ),q(θ,α,δ)]，

对于有时变时滞存在时的变分贝叶斯辨识，优点是可以将时滞作为隐变量，可得到时滞点的后验分布估计。另外，通过上述在参数先验分布中引入稀疏因子，可在VB框架下同时进行结构选择和参数估计，方便实施。

现简要说明VB参数辨识机制：

假设用C_obs，C_mis，β，Θ分别表示观测变量、缺失变量、模型结构和未知参数，我们可得到下式的log似然表达：

logp(C_obs,C_mis|β)＝logp(C_obs,C_mis,Θ|β)-logp(C_mis,Θ|C_obs,β)

向上式中加入变分后验q(C_mis,Θ)，其进行分解，q(C_mis,Θ)＝q(C_mis)q(Θ)，可得：

对上式等号两边同时对q(C_mis,Θ)取期望可得：

定义Evidence Lower Bound为：

定义KL散度为：

因此有logp(C_obs|β)＝L(q)+KL(q||p)。根据KL散度的定义，我们引入的缺失/ 隐变量和代求参数的变分后验联合分布(需要求的)q(C_mis,Θ)越接近于真实的后验分布P(C_mis,Θ|C_obs,β)，则KL(q||p)的值越小，此时L(q)越接近于 logp(C_obs|β)，使logp(C_obs|β)最大(最大似然)便等价于使L(q)最大。换句话说，我们可以通过求L(q)关于q(C_mis,Θ)的最大值，来得到q(C_mis,Θ)。求出q(C_mis,Θ)后，返代回求出的L(q)越大，则辨识的效果越好。

因为q(C_mis,Θ)可分解q(C_mis,Θ)＝q(C_mis)q(Θ)，求L(q)关于q(C_mis,Θ)的最大值，可用L(q)分别对q(C_mis)和q(Θ)求偏导并使其等于0来获得q(C_mis)和q(Θ) 的表达式。求q(C_mis)的过程叫做VB的E步；求q(Θ)的过程叫做VB的M步。下面给出q(C_mis)和q(Θ)的一般表达式：

VB参数辨识过程介绍完毕。

该步骤有两个循环，外层循环为结构修剪循环，内层循环为在特定结构下的参数循环。

外层循环初始设置:初始模型集为最大模型集

结构循环标识为s，令s＝0。

内层循环过程：(该过程在结构

下进行，为方便起见，省略s标识)

1)收集辨识输入输出数据，给出超参数{a₀,b₀,c₀,d₀,D}初始值,按照先验分布初始化参数

其中

是以

为对角线值的对角矩阵，

符号“-”代表平均值，上标表示参数循环标识，此处为0，表示初始值，M是当前模型集的子模型个数，选择一个小的正数ε，作为参数收敛标准；

2)使用VB方法进行参数辨识，首先用VB的E步得到时滞λ的后验分布更新公式

其中，

表示在已知Φ_k，

λ_k＝j下，第k个时刻输出为y_k的概率；

表示在第t次迭代时，时滞为j的概率；y_k表示k时刻的输出值，Φ_kj表示k时刻时，时滞为j时的模型集；

为t次参数迭代下θ的均值；

表示θθ^T在θ的后验分布下的值， 3)在上步中得到的q(λ_k＝j)下进行VB的M步,得到参数

Var(θ)^t+1，θ^t+1，

和

的更新公式

按下列公式得到参数

Var(θ)^t+1和θ^t+1的更新公式：

按下列公式得到参数

的更新公式：

q(δ)＝Gamma(c,d)

按下列公式得到

的更新公式

其中，

是

中的第m个值，

是Var(θ)^t中的第m行m列的值；

按下列公式得到参数

的更新公式：

4)按下列公式判断参数是否收敛

若满足上式条件，则跳出内层循环得到参数

Var(θ)^t+1、θ^t+1、a 、

c、d^t+1、q(λ_k＝j)^t+1和

若不满足上式，则t＝t+1返回步骤2)

5)参数收敛后，在当前结构和参数下计算Evidence Lower Bound, F[q(λ),q(θ,α,δ)]的计算公式为：

F[q(λ),q(θ,α,δ)]＝f₁-f₂+f₃-f₄，

其中，

外层循环过程

1)修剪模型集

计算修剪子模型的指标ARD^s，由内层循环得到的参数

取逆表示，记为

修剪条件：log(ARD^s)中的值小于限值

的冗余子模型会被从模型中修剪出去，其中，限值

的定义如下所示：

修剪后得到新的模型集结构

模型集的大小M也得到了更新；

2)判断模型个数是否为1，若M＝1，则退出外层循环；否则令s＝s+1，并

在当前修剪后的模型结构下返回内层循环初始步骤，再次进行参数估计。步骤三：在步骤二的基础上，在最终只保留了一项子模型的结构时，寻找不同结构参数对应的F[q(λ),q(θ,α,δ)]^s(s＝0,1,2,…是模型修剪的次数)中的最大值，此时最大的F[q(λ),q(θ,α,δ)]函数值对应的模型结构和参数即为最优模型，完成结构辨识和参数辨识。

在最终只保留了一项子模型的结构时，寻找不同结构参数对应的F[q(λ),q(θ,α,δ)]^s(s＝0,1,2,…是结构循环标识)中的最大值，可见：

此时最大的F[q(λ),q(θ,α,δ)]^s函数值对应的模型结构

和参数Θ*即为最优，完成结构辨识和参数辨识。

有益效果：

由于不确定实验室分析或网络拥堵，时滞在每次采样时变化。无需知道时滞的精确范围，只需给定时滞上界即可。在VB辨识下，可得到每个采样时刻时滞的分布估计而不是点估计。

结构和参数的辨识可在VB框架下迭代进行，在每次模型结构修剪过程中，通过在参数先验中引入稀疏因子，实现了每个子模型与数据预测的相关性度量。对于给定的数据集，不同子模型的显著性是可以直接比较的，因此提供了一种快速简单的模型选择方法。

仿真结果表明，该算法对参数初值不敏感，参数值设置为0也可收敛到真实值。

附图说明

图1为本发明辨识流程图。

图2为一个水箱模型

图3为模型结构修剪图

图4为时滞估计图

具体实施方式

参照附图并举实施例说明本发明的具体实施方式。

图1为本发明提供的一种针对时变时滞下NARX模型结构和参数辨识流程图。

以一个水箱模型为例，见图2，u_k和

分别为阀的开度和进水流量。u_k和

之间存在的非线性关系为：

y_k为水箱2的液位高度。以

为输入、y_k为输出的线性部分的传递函数为： G(z)＝(b₁z^-1+b₂z^-2)/(1+a₁z^-1+a₂z^-2)。系统的真实参数为 [a₁,a₂,b₁,b₂,ζ₁,ζ₂]^T＝[-0.5 0.83 0.36 1.1 1 0.58]^T。

输入序列{u_k}取为不相关的持久激励序列，Δt＝10s更新一次。由压力传感器得到的输出序列{y_k}以相同的时间间隔更新。我们对输入信号手动施加一个变化的时间延迟为10s，20s和30s，以模拟由测量设备和信号通信引入的可能的时间延迟。因此仿真时，实际的时滞为{1,2,3}，并且三种时滞时间按概率β₁＝0.2, β₂＝0.5，β₂＝0.3产生。{v_k}是方差为0.04的高斯白噪声序列。取N＝600个数据进行仿真。

第一阶段：为方便说明，设置n_u＝n_y＝2，n_l＝3，因此得到个子模型个数为 M＝34的初始模型集。建立NARX的多项式子模型加权和形式为：

第二阶段：为进行结构辨识引入稀疏因子α。设置参数先验，将θ的先验分布设置为正态分布；α和δ的分布设置为Gamma分布。

p(θ|α)＝p(θ|0,(A)^-1)

p(δ)＝Gamma(c₀,d₀)

假设关于时滞先验信息为最大时滞时间不超过D，因此设时滞j∝{1,2,…,D}

超参数初始化为：{a₀,b₀,c₀,d₀,D}＝{0.01,0.0001,0.01,0.0001,4}

初始结构下的参数初始化为：

参数辨识收敛后的子模型对应的权重值以及是否修剪子模型项如下表所示：

计算当前结构和参数下的Evidence Lower Bound，计算出的修剪限值为：

第一轮结构修剪后，新的子模型集为：

使用新的模型集重新进行参数辨识，参数收敛后再次计算Evidence Lower Bound和限值

第二次模型修剪过程见下表：

第二次模型修剪后得到的结构为:

用新模型继续进行参数辨识，直至模型修剪到剩一个。此时s＝6

在每种模型下计算得到的Evidence Lower Bound为

s	Evidence Lower Bound
		0	-INF
1	517.513701
		2	536.349211
3	-524.065867
		4	-3139.817458
5	-3535.626748
		6	-3653.818099

当s＝2时，Evidence Lower Bound有最大值。此时模型的结构和参数为：

最优模型结构下，时滞参数的估计值见下表：

结构辨识结果如图3所示，可知第三次结构迭代，s＝2时即可将正确模型选择出来。

时滞辨识结果如图4所示，时滞正确率为：0.825。

Claims (4)

1.一种在时变时滞下NARX模型结构和参数的变分贝叶斯辨识方法，其特征在于包括两大阶段，

第一阶段，对水箱模型进行结构预设，将NARX非线性模型转化为多项式子模型加权和的形式；

第二阶段，对多项式子模型加权和的形式进行结构辨识和参数辨识，即选出能够表达系统的最佳多项式子模型集Φ_k，以及求出时变时滞下模型对应的参数{Θ,β}，其中Θ＝{θ,δ,α}，β＝{β₁,β₂,…,β_j,…}，θ表示子模型的权重向量，δ表示噪声方差的逆，α表示参数θ的先验分布中引入的稀疏因子，β_j表示时滞值j出现的概率，具体包括以下步骤：

步骤一：设置参数的先验分布；

步骤二：在VB框架下，在步骤一的基础上进行特定结构下的参数辨识，再用其中的稀疏因子对多项式子模型加权和的形式进行结构修剪即把冗余的多项式子模型剔除，用修剪后的新结构再次进行参数辨识，直到多项式子模型加权和的形式最终只保留一项子模型；计算每种结构和及其对应参数下的Evidence Lower Bound，即F[q(λ),q(θ,α,δ)]，

所述的步骤二为：

该步骤有两个循环，外层循环为结构修剪循环，内层循环为在特定结构下的参数循环，

外层循环初始设置:初始模型集为最大模型集

结构循环标识为s，令s＝0.

内层循环过程：该过程在结构

下进行，为方便起见，省略s标识；

1)收集辨识输入输出数据，给出超参数{a₀,b₀,c₀,d₀,D}初始值,按照先验分布初始化参数

其中

是以

为对角线值的对角矩阵，

符号“-”代表平均值，上标表示参数循环标识，此处为0，表示初始值，M是当前模型集的子模型个数，选择一个小的正数ε，作为参数收敛标准；

2)使用VB方法进行参数辨识，首先用VB的E步得到时滞λ的后验分布更新公式

其中，

表示在已知Φ_k，

λ_k＝j下，第k个时刻输出为y_k的概率；

表示在第t次迭代时，时滞为j的概率；y_k表示k时刻的输出值，Φ_kj表示k时刻时，时滞为j时的模型集；

为t次参数迭代下θ的均值；

表示θθ^T在θ的后验分布下的值，

3)在上步中得到的q(λ_k＝j)下进行VB的M步,按下列公式得到参数

Var(θ)^t+1，θ^t+1，

和

的更新公式

按下列公式得到参数

的更新公式：

按下列公式得到

的更新公式

其中，

是

中的第m个值，

是Var(θ)^t中的第m行m列的值；

按下列公式得到参数

的更新公式：

4)按下列公式判断参数是否收敛

若满足上式条件，则跳出内层循环得到参数

Var(θ)^t+1、θ^t+1、a、

c、d^t+1、q(λ_k＝j)^t+1和

若不满足上式，则t＝t+1返回步骤2)

5)参数收敛后，在当前结构和参数下计算Evidence Lower Bound,F[q(λ),q(θ,α,δ)]的计算公式为：

F[q(λ),q(θ,α,δ)]＝f₁-f₂+f₃-f₄，

其中，

外层循环过程

1)修剪模型集

计算修剪子模型的指标ARD^s，由内层循环得到的参数

取逆表示，记为

修剪条件：log(ARD^s)中的值小于限值

的冗余子模型会被从模型中修剪出去，其中，限值

的定义如下所示：

修剪后得到新的模型集结构

模型集的大小M也得到了更新；

2)判断模型个数是否为1，若M＝1，则退出外层循环；否则令s＝s+1，并在当前修剪后的模型结构下返回内层循环初始步骤，再次进行参数估计；

步骤三：在步骤二的基础上，在最终只保留了一项子模型的结构时，寻找不同结构参数对应的F[q(λ),q(θ,α,δ)]^s中的最大值，此时最大的F[q(λ),q(θ,α,δ)]函数值对应的模型结构和参数即为最优模型，完成结构辨识和参数辨识，s＝0,1,2,…是模型修剪的次数。

2.根据权利要求1所述的一种在时变时滞下NARX模型结构和参数的变分贝叶斯辨识方法，其特征在于，第一阶段所述的带时变时滞的单输入单输出NARX模型描述如下：

其中，f(·)是某个非线性函数；{u_k}和{y_k}分别是采样得到的输入和输出数据，即分别是阀的开度和水箱的液位高度；下标k表示第k个采样时刻；n_u,n_y分别是输入和输出的最大动态阶次；λ_k为第k个采样时刻出现的时滞，

{v_k}是均值为0，方差为δ^-1的高斯分布白噪声；

可对动态阶次n_u、n_y及多项式阶次n_l设置上限，将NARX模型表示为多项式子模型加权和的形式：

Φ_k＝[φ₁,φ₂,…,φ_M]，θ＝[θ₁,θ₂,…,θ_M]^T

其中，φ_m是由

构成的多项式子模型，M表示多项式子模型的个数，θ_m表示每个子模型的权重参数，因此得一个初始的模型集Φ_k。

3.根据权利要求1所述的一种在时变时滞下NARX模型结构和参数的变分贝叶斯辨识方法，其特征在于，所述的步骤一为：

所述的设置参数的先验分布指对使用变分贝叶斯方法进行辨识时所需要的参数设置参数先验，包括设置参数θ的先验分布，α的先验分布，噪声方差的逆δ的先验分布，以及时滞的先验分布，具体如下：

为实现结构辨识，在参数θ的先验分布中引入稀疏因子α，同时为实现共轭分布，选择θ的先验分布为正态分布，具体如下，

p(θ|α)＝p(θ|0,(A)^-1)

其中，α＝(α₁,α₂,…,α_m,…,α_M)^T＝diag(A)，α_m用于控制θ_m变化的幅度，由分析可知，如果

则θ_m越接近于0，说明第m个子模型项对数据来说不显著，为冗余结构，应该将其从模型集中移除；

同样，为得到一个共轭分布，将α的先验分布设置为Gamma分布：

考虑到辨识问题的一般情况，噪声的方差未知，同样将方差的逆δ作为需要辨识的参数并设置先验分布，即如下：

p(δ)＝Gamma(c₀,d₀)

假设时滞出现的最大值为D，设置时滞的先验分布，将每种时滞出现的概率初始化为相同的值β_j即如下：

其中，{a₀,b₀,c₀,d₀,D}是需要人为给定的初始化超参数。

4.根据权利要求1所述的一种在时变时滞下NARX模型结构和参数的变分贝叶斯辨识方法，其特征在于所述的步骤三为：

在最终只保留了一项子模型的结构时，寻找不同结构参数对应的F[q(λ),q(θ,α,δ)]^s中的最大值，s＝0,1,2,…是结构循环标识，可见：

此时最大的F[q(λ),q(θ,α,δ)]^s函数值对应的模型结构

和参数Θ^*即为最优，完成结构辨识和参数辨识。

Images