基于变分贝叶斯算法的线性变参数系统辨识

基于变分贝叶斯算法的线性变参数系统辨识
李寒霜;赵忠盖;刘飞
【摘要】线性变参数系统(LPV)将多阶段、非线性的过程建模转化为线性多模型的辨识问题,是解决非线性过程建模的一个有效手段.由于实际工业过程存在各种干扰因素,导致被建模系统呈现随机性及模型参数的不确定性.针对这一问题,考虑采用变分贝叶斯(VB)算法对LPV模型进行辨识.该算法首先给定参数相应的先验分布,通过最大化目标函数的下界,从而估计得到参数的后验分布.不仅可实现对参数的点估计,同时量化了估计值的不确定性.针对典型二阶过程和连续搅拌反应釜(CSTR),运用提出的算法进行仿真实验,表明了该贝叶斯估计方法的优越性.
【期刊名称】《化工学报》
【年(卷),期】2018(069)007
【总页数】10页(P3125-3134)
【关键词】非线性过程;线性变参数系统;多模型;变分贝叶斯算法;参数估计
【作者】李寒霜;赵忠盖;刘飞
【作者单位】江南大学轻工过程先进控制教育部重点实验室,江苏无锡 214122;江南大学轻工过程先进控制教育部重点实验室,江苏无锡 214122;江南大学轻工过程先进控制教育部重点实验室,江苏无锡 214122
【正文语种】中文
【中图分类】TP301.6
引言
实际工业过程大多呈现显著的非线性特性，并且机理复杂，工况多变，很难采用单一模型描述过程的全局动态行为。

线性变参数 (linear parameter varying，LPV) 模型基于简单的线性结构，通过参数的变化即可精确描述复杂的非线性时变过程，得到了广泛关注[1-6]。

LPV方法最初由Shamma等[7]提出，用于增益调度控制器设计的模型，参数为可测或可计算时变参数的函数。

该时变信号通常被称为调度变量，反映了过程的工作模态[8]。

近年来，关于LPV模型的辨识方法逐渐成为研究焦点。

Xu等[9]考虑在系统典型工作点处建立局部模型，再通过权函数将各个子模型融合，并提出可以采用预测误差模型辨识方法、子空间辨识方法等得到局部线性模型。

但是，该方法采用三次样条函数作为权函数，而确定三次样条函数合适的阶数具有一定难度；Jin 等[10]将期望最大化算法 (EM) 引入LPV模型的辨识中，通过极大似然的方法辨识模型中未知参数，得到了较高的辨识精度；在Jin等的研究基础上，Yang等[11]针对含有量测时滞的非线性系统，采用 EM 及其广义算法辨识得到了系统的LPV 模型。

而在实际工业过程中，由于恶劣环境中存在各种干扰因素，所有驱动系统行为的变量难以完全观测，从而导致被建模系统呈现随机性及模型参数的不确定性[12]。

但是，上述辨识方法均没有考虑到这一问题，只是对参数进行单点估计，导致模型精度的降低。

为解决上述问题，近年来，一种能适应高维观测环境，并将未知参数视为随机变量的参数估计方法逐渐被人们接受。

Attias[13]最早提出将变分近似的思想应用于图像模型的参数推断中，成为了VB方法的雏形；Yang等[14]针对存在缺失数据的线性FIR模型，基于VB算法推导得到参数的贝叶斯估计；Zhao等[15]研究了存在异常点数据且含有不确定量测时滞的线性ARX模型，基于VB算法提出了模型
的鲁棒辨识问题；此外，VB算法已被应用在其他研究领域，如狄利克雷混合模型[16]和连续离散随机动态系统[17]，且都取得了不错的成果。

文献[15]研究了线性系统的参数辨识问题，即采用VB算法实现对ARX模型的参数估计。

考虑到实际工业过程大多为非线性系统，本文提出一种基于VB算法的LPV模型辨识方法，不仅需要对多个线性子模型进行参数辨识，还需对加权函数进行估计，比单个线性模型的辨识过程更为复杂。

该方法采用一个形式较为简单的自由分布描述 LPV参数的分布，通过解析近似的推理方法，不断最大化边缘似然函数的下界，迭代地更新变分参数，直至近似分布足够逼近参数真实的后验分布，从而实现对LPV模型参数的估计。

最后，将方法应用于仿真算例中，以证明其有效性。

1 线性变参数模型
LPV模型建立的基本方法是应用简单的局部模型来表示系统在不同工作情况下的特性，然后通过权函数将若干局部模型融合为全局模型，即LPV模型是由一组线性模型通过权函数加权融合所得。

由于 ARX模型具有较好地逼近线性动态系统的性能[18]，因此本文采用 ARX模型作为局部线性模型结构，表示如下
式中，变量Ik代表k时刻变量隶属的局部模型，取值为1,2,…,i,…,M 中的一个。

xk代表系统在第 k次采样时的回归量，ek是均值为0、方差为的高斯分布噪声，未知。

θIk表示局部模型的参数向量，同时引入参数ηi表示局部模型参数的精度[19]。

yk和uk分别为系统的输出和输入，na和nb分别表示输出和输入的阶数。

不失一般性，本文取 SISO系统，且输出输入阶数na=nb=2。

假定系统中有M个工作点，即选择M个局部模型来描述系统的动态特性。

为使各子模型平滑融合，采用能兼顾模型精度和辨识可靠性的高斯分布函数为权重函数[20-22]，表示如下
其中，k=1,2,…,N，i=1,2,…,M。

H1:N表示系统的调度变量，表示系统的主要工作点，Oi表示第i个局部模型的有效宽度。

归一化后的权值函数代表在k时刻第Ik=i个局部模型起作用的概率
那么，对于系统输出的预测如式(5)所示
式中，表示k时刻LPV模型的预测输出，表示k时刻第i个ARX子模型的预测输出。

在整个过程中需要估计的参数包含子模型参数θi，子模型参数的精度ηi，模型噪声的精度σi以及子模型的有效宽度Oi。

因此，LPV模型辨识问题就转化为基于测得的辨识数据采用VB算法估计模型中的未知参数。

2 VB算法
VB算法针对模型中的隐藏和未知变量，通过最大化变分参数的对数边缘似然函数的下界(目标函数)，可以得到隐藏变量和未知参数的贝叶斯估计。

具体实现方法如下。

令系统观测数据集为Cobs，隐藏数据集为Cmis，待估计的参数为Θ。

在贝叶斯估计中，观测数据的边缘似然函数可由式(6)计算
式(6)中有难以处理的高维复杂积分，VB算法通过引入联合分布q(Cmis,Θ)，也称变分后验分布，来近似计算隐藏变量和未知参数的后验分布p(Cmis,Θ)。

进而，对数边缘似然函数为
其中，在式(7)的推导中，应用了Jensen不等式。

在 VB的框架下，假定联合分布
q(Cmis,Θ)是可分解的，即q(Cmis,Θ)=q(Cmis)q(Θ)。

因此，式(7)表示为
将对数边缘似然函数 ln p(Cobs)和下界函数F[q(Cmis),q(Θ)]做差，可得
其中，KL(q||p)是 Kullback-Leibler(KL)散度[23-24]，表示联合分布q(Cmis,Θ)与实际分布p(Cmis,Θ|Cobs)之间的差异。

由式(9)可知，最小化 KL(q||p)等价于最大化下界函数F[q(Cmis),q(Θ)]。

因此，VB算法通过最大化下界函数，迭代地更新隐变量和未知参数的后验分布表达式，直至算法收敛，即可得到真实后验分布p(Cmis,Θ|Cobs)的近似解
q(Cmis,Θ)。

3 基于VB算法的参数估计
3.1 辨识方法总体流程
在第1部分建立的LPV模型中，观测到的数据为缺失数据为 Cmis={I1:N}，需要
估计的参数集合为{θ1:M, η1:M, σ1:M,O1:M}。

在VB算法中，通过选择合适的先验分布来表征参数Θ={θ1:M, η1:M, σ1:M}的不确定性，从而求得它们的后验分布。

而对于局部模型的有效宽度O，是直接求取它的点估计，因此将参数Θ的求解与
O分离开。

基于VB算法对LPV模型的辨识方法可归纳为以下步骤。

（1）初始化：给模型中未知参数Θ分配合适的先验分布，并初始化参数向量Θ、O及先验分布中的超参数；
（2）E步：固定参数集，关于隐变量最大化下界函数，得到隐变量I后验分布的
更新式q(I)；
（3）M步：固定隐变量，同时考虑到各参数的独立性，即q(Θ)=q(θ)q(η)q(σ)，则关于每个参数分别最大化下界函数，得到各自后验分布的更新式q(θ)、q(η)和
q(σ)，并计算得到O的点估计；
（4）估计下界函数F[q(I),q(Θ)]；
（5）迭代过程：将获得新的估计值代入步骤(2)，不断重复步骤(2)～步骤(4)，直
至下界函数达到最大值。

下面对各步骤进行详细介绍。

3.2 参数的先验分布
模型参数Θ的联合先验分布可以表示为
采用如下的共轭先验分布[15]，共轭先验分布的使用确保了参数后验分布与先验分布的一致性，使得贝叶斯分析得到了极大的简化
其中，表示(na+nb)维的单位矩阵，下文均简记为D；a0、b0、c0、d0是先验Gamma分布的超参数。

3.3 E步
在E步中，固定q(Θ)，关于q(I)最大化下界函数，下界函数F[q(I),q(Θ)]表示如下
通过求解如下的优化问题
可以得到q(I)的表达式
其中〈·〉q(Θ)表示关于Θ 的期望运算。

式(16)、式(17)的详细推导过程见附录1。

根据概率链式法则，全部数据的对数似然函数可以扩展为
其中，为常数。

因此，式(17)可以表示成
即那么式(16)中的q(I)可表示如下
同时，考虑到可以得到q(Ik)
为简化表达式，令
则
接下来，需要求解出Aki的表达式。

式(18)中，表示的是ARX模型分布概率，符合高斯分布，根据高斯分布概率函数可以得到
表示权重函数，可以根据式(4)求得。

式(18)可表示如下
进一步，可得到式(27)中的期望对数似然
将推导得到的式(27)代入式(22)，可得到隐藏变量Ik的概率分布q(Ik)。

3.4 M步
在M步中，固定q(I)，关于q(Θ)最大化下界函数F[q(I),q(Θ)]。

假定q(Θ)可分解为如下的形式，可以对分解后的每一项分别最大化下界函数
在该步中，下界函数可表示为
其中表示如下
首先关于q(θi)最大化下界函数，得到如下的优化问题
为求解式(32)，需计算关于q(θi)拉格朗日函数的一阶导数
其中λ是拉格朗日乘子。

从而，可以得到q(θi)为
由式(34)推导可得，q(θi)是高斯分布的密度函数，即q(θi)=N(θi|E(θi),Var(θi))，且该高斯分布的均值和方差表示如下
下一步采用同样的方法关于q(ηi)最大化下界函数，可以得到
其中Cη是与ηi无关的常数。

由式(37)可得，q(ηi)是Gamma分布的密度函数，即且
根据Gamma分布的性质可得到ηi的期望值为
同样地，可以得到q(σi)的表达式也是 Gamma分布的密度函数，即且
σi的期望值可计算如下
对于局部模型有效宽度，应用寻优的方法来搜索最优解，其数学表达式为
本文通过一种非线性数值优化方法来计算Oi，即直接采用 Matlab中有约束的非
线性优化函数fmincon进行求解[25]。

3.5 估计下界函数
运用每次迭代后更新的估计值 q(Ik)、q(θi)、q(ηi)、q(σi)和 Oi计算下界函数，观察是否达到收敛条件。

下界函数F[q(I),q(Θ)]的表达式见式(14)，直接计算式中的积分比较困难，考虑通过其他数学运算将积分消去[26]。

式(14)的第1项可以转换为
式(45)的具体推导见附录2。

式(14)的第2项可表示为
因此，下界函数F可表示如下
4 仿真验证
4.1 数值仿真
以模型参数变化的二阶过程为例，验证前述LPV模型辨识方法，二阶过程的传递函数如式(48)所示
其中，过程增益 K(h)和过程时间常数τ(h)均为调度变量h的非线性函数，具体表示如下
可以看出，在整个过程的操作范围内，过程增益和时间常数的变化比较剧烈，单一的线性模型难以描述整个工作范围内的动态特性，因此采用LPV的方法对该过程建模。

选取工作点为整个过程中调度变量的轨迹及实际的输入输出数据如图1所示，图2表示归一化后的权重与调度变量之间的关系。

将本文提出的辨识方法应用于二阶非
线性过程，在不同操作条件下辨识局部模型，下界函数值F随迭代次数的变化曲线如图3所示，显然算法在7次迭代内已达到收敛。

图1 二阶非线性过程训练数据集Fig.1 Training data set of second-order nonlinear process
图2 二阶非线性过程不同操作点下的局部模型权重Fig.2 Normalized weight of each local model
图3 下界函数值F与迭代次数的变化曲线Fig.3 Calculated lower bound versus iteration times
为进一步说明VB和EM算法的区别，采用文献[10]中提出的基于EM算法的LPV 模型辨识方法对该过程仿真。

表1列出了ARX子模型的真实参数值以及分别基于VB和EM算法估计得到的模型参数值，可以看出应用本文提出算法估计得到的参数结果更接近真实值。

同时，以子模型1为例，基于VB算法估计得到的参数概率分布函数曲线如图4所示。

图5为模型自我验证的结果，为方便比较，图中仅给出了区间[400,500]内的数据。

同时，模型的质量通过相关误差量测得[27]
表1 分别基于VB算法和EM算法的参数估计结果Table 1 Estimated parameter results based on VB and EM methodLocal model Parameter True value Estimated value based on VB method Estimated variance based on VB method Estimated value based on EM method 1 a1 1.5088 1.5565
1.0817×10−4 1.6282 a2 −0.7515 −0.7745 1.0217×10−4 −0.8016 b1 0.5062 0.5222 1.3996×10−4 0.5435 b2 −0.2635 −0.3128 1.4664×10−4 −0.4001 2 a3 1.8240 1.8460 0.5813×10−4 1.8713 a4 −0.9117 −0.9186 0.5562×10−4 −0.9363 b3 0.6403 0.6526 1.0777×10−4 0.6218 b4 −0.5527 −0.5990
1.1204×10−4 −0.5763 3 a5 1.9437 1.9442 0.7387×10−4 1.9235 a6 −0.9718
−0.9761 0.7339×10−4 −0.9603 b5 0.4795 0.5441 1.5500×10−4 0.4858 b6 −0.4514 −0.5073 1.5669×10−4 −0.4568
其中Var定义为信号的方差，y为真实输出，为模型输出。

VB算法和EM算法自我验证的相关误差分别为 0.19%和 0.34%。

为更好地验证模型的有效性，选取另一组输入输出数据进行交叉验证，VB和EM算法交叉验证的结果如图6所示，计算得到的相关误差分别为0.46%和1.03%。

结果表明，与EM算法相比，采用VB 算法对模型参数的估计精度更高。

图4 局部模型1的参数分布估计Fig.4 Parameter distribution estimations of local model 1
图5 二阶非线性过程自我验证结果Fig.5 Self-validation result of second-order nonlinear process
图6 二阶非线性过程交叉验证输出结果Fig.6 Cross-validation result of second-order nonlinear process
4.2 CSTR仿真实验
CSTR是化工、发酵、石油生产等工业生产过程中应用最广泛的一种化学反应器。

不可逆放热反应发生在体积不变的反应器中，通过冷却剂冷却反应器中的温度。

图7为一个装有冷却管的CSTR装置示意图。

该过程的第一原理模型如式(52)所示，系统模型参数及稳态值见表2[28-30]。

图7 CSTR系统装置示意图Fig.7 Diagrammatic sketch of CST R
表2 CSTR模型参数及稳态值Table 2 CSTR model parameters and their steady state valuesNote: 1 cal=4.1868 J.Parameter Value production concentration of component A,CA output 1 temperature of reactor,T output 2 coolant flow rate, qc input feed temperature, T0/K 350.0 specific heats, cp, cpc/(cal·g−1·K−1) 1 liquid density, ρ, ρc/(g·L−1) 1×103 heat of
reaction, −ΔH/(cal·mol−1) −2×105 reaction rate constant, k0/min−1
7.2×1010 activation energy term, (E/R)/K 1×104 heat transfer term,
hA/(cal·min−1·K−1) 7×105 inlet coolant temperature, Tc0/K 350 process flow rate, q/(L·min−1) 100 reactor volume, V/L 100 feed concentration of component A, CA0/(mol·L−1) 1
本文中，反应器内生产物浓度CA为主要输出，冷剂流量qc为输入变量。

由于冷剂流量对过程动态有重要影响，因此选择它为调度变量，范围是97～105 L·min−1，且预先选定 97、100、101、103、105 L·min−1 5个工作点[27]。

输入信号为幅值在区间[97,105]内的随机二进制信号(RBS)，在获取的输出数据中加入方差为0.015的过程白噪声，系统的输入输出数据以及整个操作范围内调度变量的值如图8所示。

图9为归一化后的权重曲线。

将本文提出的辨识方法和EM 方法分别应用于CSTR系统的仿真数据中，自我验证的相关误差分别为 3.14%和3.43%。

进一步，选取另一组输入输出数据构成交叉验证集来检验辨识模型的预测性能，计算得到 VB和 EM 方法交叉验证的相关误差分别为 3.28%和3.52%。

自我验证和交叉验证的结果如图10和图11所示，从图中可以看出本文所述方法建立的模型拟合精度较高，从而验证了辨识方法的有效性。

图8 CSTR系统训练数据集Fig.8 Training data set of CSTR system
图9 CSTR系统不同操作点下的局部模型权重Fig.9 Normalized weight of each local model
图10 CSTR系统自我验证输出结果Fig.10 Self-validation result of CSTR system
图11 CSTR系统交叉验证输出结果Fig.11 Cross-validation result of CSTR system
5 结论
本文提出了一种基于VB算法的LPV模型辨识方法。

该算法以概率分布的形式描
述参数的不确定性，将模型中多变量的联合概率分布近似为各变量边缘概率分布的乘积，可迭代估计得到各参数的后验分布，该后验分布能够将估计参数的统计特性，包括均值和方差解析地表示出来。

与忽略参数不确定性的EM算法相比，VB算法在参数估计方面展示了更优异的性能。

最后对二阶非线性过程和CSTR系统分别进行仿真验证。

结果表明，该方法可以有效、精确地实现对LPV模型参数的估计。

References
附录1 式(16)、式(17)推导过程
由式(14)可得
计算关于q(I)拉格朗日函数的一阶导数
从而
将式(A4)代入约束中, 可得
为简化q(I)的表达式, 令
可以得到q(I)的更新表达式
式(16)、式(17)即可推导得出。

附录2 式(45)推导过程
式(45)即可推导得出。

【相关文献】
[1]VERDULT V, VERHAEGEN M. Subspace identification of multivariable linear parameter-varying systems[J]. Automatica, 2002,38(5): 805-814.
[2]TOTH R, HEUBERGER P S C, HOF P M J V D. Asymptotically optimal orthonormal basis functions for LPV system identification [J].Automatica, 2009, 45(6): 1359-1370.
[3]ZHAO Y, HUANG B, SU H Y, et al. Prediction error method for identification of LPV models[J]. Journal of Process Control, 2012,22(1): 180-193.
[4]TEHRANI E S, KEARNEY R E. A non-parametric linear parameter varying approach for identification of linear time-varying systems[J].IFAC-Papers Online, 2015, 48(28): 733-738.
[5]HOFFMANN C, WERNER H. A survey of linear parameter-varying control applications validated by experiments or high-fidelity simulations[J]. IEEE Transactions on Control Systems Technology,2015, 23(2): 416-433.
[6]ZHAO P, NAGAMUNE R. Switching LPV controller design under uncertain scheduling parameters[J]. Automatica, 2017, 76: 243-250.
[7]SHAMMA J S, ATHANS M. Guaranteed properties of gain scheduled control for linear parameter-varying plants[J]. Automatica, 1991,27(3): 559-564.
[8]ZHU Y C, XU Z H. A method of LPV model identification for control[J]. IFAC Proceedings Volumes, 2008, 17(1): 5018-5023.
[9]XU Z H, ZHAO J, QIAN J X, et al. Nonlinear MPC using an identified LPV model[J]. Industrial & Engineering Chemistry Research, 2009,48(6): 3043-3051.
[10]JIN X, HUANG B, SHOOK D S. Multiple model LPV approach to nonlinear process identification with EM algorithm[J]. Journal of Process Control, 2011, 21(1): 182-193. [11]YANG X Q, HUANG B, GAO H J. A direct maximum likelihood optimization approach
to identification of LPV time-delay systems[J].Journal of the Franklin Institute, 2016, 353(8): 1862-1881.
[12]GOODFELLOW L, BENGIO Y, COURVILLE A. Deep Learning[M].Massachusetts: The MIT Press, 2016: 54-57.
[13]ATTIAS H. A variational Bayesian framework for graphical models[J].International Conference on Neural Information Processing Systems,1999, 12: 209-215.
[14]YANG X Q, YIN S. Variational Bayesian inference for FIR models with randomly missing measurements[J]. IEEE Transactions on Industrial Electronics, 2016, 64(5): 1-9.
[15]ZHAO Y J, FATEHI A, HUANG B. Robust estimation of ARX models with time varying time delays using variational Bayesian approach[J].IEEE Transactions on Cybernetics, 2017, 99: 1-11.
[16]MA Z, RANA P K, TAGHIA J, et al. Bayesian estimation of Dirichlet mixture model with variational inference[J]. Pattern Recognition,2014, 47(9): 3143-3157.
[17]ALA-LUHTALA J, SARKKA S, PICHE R. Gaussian filtering and variational approximations for Bayesian smoothing in continuousdiscrete stochastic dynamic systems[J]. Signal Processing, 2015, 111:124-136.
[18]BACHNAS A A, TOTH R, LUDLAGE J H A, et al. A review on datadriven linear parameter-varying modeling approaches: a high-purity distillation column case study[J]. Journal of Process Control, 2014,24(4): 272-285.
[19]LU Y J, KHATIBISEPEHR S, HUANG B. A variational Bayesian approach to identification of switched ARX models[J]. IEEE Conference on Decision and Control, 2014, 2015: 2542-2547.
[20]HUANG J Y, JI G L, ZHU Y C, et al. Identification of multi-model LPV models with two scheduling variables[J]. Journal of Process Control, 2012, 22(7): 1198-1208.
[21]YOU J, LU J G, ZHU Y C, et al. Identification of multimodel LPV models with asymmetric Gaussian weighting function[J]. Journal of Applied Mathematics, 2013,
2013(2): 1-12.
[22]YOU J, YANG Q M, LU J G, et al. Identification of LPV models with non-uniformly spaced operating points by using asymmetric Gaussian weights[J]. Chinese Journal of Chemical Engineering, 2014, 22(6):795-798.
[23]BERNARDO J M, BAYARRI M J, BERGER J O, et al. The variational Bayesian EM algorithm for incomplete data: with application to scoring graphical model structures[J]. Bayesian Statistics, 2003,78(78): 417-429.
[24]LU Y J, HUANG B, KHATIBISEPEHR S. A variational Bayesian approach to robust identification of switched ARX Models[J]. IEEE Transactions on Cybernetics, 2017, 46(12): 3195-3208.
[25]LU Y J, HUANG B. Robust multiple-model LPV approach to nonlinear process identification using mixture t distributions[J].Journal of Process Control, 2014, 24(9): 1472-1488.
[26]TAKEKAWA T, FUKAI T. A novel view of the variational Bayesian clustering[J]. Neurocomputing, 2009, 72(13/14/15): 3366-3369.
[27]CHEN L, TULSYAN A, HUANG B, et al. Multiple model approach to nonlinear system identification with an uncertain scheduling variable using EM algorithm[J]. Journal of Process Control, 2013, 23(10):1480-1496.
[28]DENG J, HUANG B. Identification of nonlinear parameter varying systems with missing
output data[J]. AIChE Journal, 2012, 58(11):3454-3467.
[29]XIONG W L, YANG X Q, HUANG B, et al. Multiple-model based linear parameter varying time-delay system identification with missing output data using an expectation-maximization algorithm[J]. Industrial& Engineering Chemistry Research, 2014, 53(27): 11074-11083.
[30]YANG X Q, YIN S. Robust global identification and output estimation for LPV dual-rate systems subjected to random output timedelays[J]. IEEE Transactions on Industrial Informatics, 2017,13(6): 2876-2885.。