基于随机分布理论的递推子空间辨识

于 淼 刘建昌 郭 戈

(1.东北大学秦皇岛分校控制工程学院,河北秦皇岛 066004;2.东北大学信息科学与工程学院,辽宁沈阳 110819;3.东北大学流程工业综合自动化国家重点实验室,辽宁 沈阳 110819)

1 引言

近年来,子空间辨识作为系统辨识的重要分支在控制领域受到了广泛的关注[1–3].针对离线模型的辨识问题,子空间辨识方法已具有较好的性能.与离散时间系统的辨识相比较,连续时间系统的辨识在工业工程、航空航天、生物学、医学及社会经济学等众多领域应用广泛[4].通过连续时间系统模型的辨识,可以较准确地判断复杂的连续时间系统的现象和过程的内在规律.此外,多数物理现象都具有连续属性,描述它们的数学模型是微分方程[5].因此,针对连续时间系统的辨识问题进行深入研究,具有重大的理论意义和实际价值.

采用子空间辨识方法建立的离线模型不能有效跟踪系统的动态变化.利用输入输出数据的更新不断地更新所建立的模型,这样才能及时地递推更新模型.递推辨识可以不断地从系统获得更新数据,从而提高辨识的模型精度直至达到辨识要求.递推辨识过程中主要使用正交三角分解(QR decomposition,QRD)和奇异值分解(singular value decomposition,SVD)等线性代数工具,虽然提高了算法的数值鲁棒性,但也相应增加了子空间辨识的在线递推困难[6–7].所以如何降低递推子空间辨识方法的计算量是首要问题.文献[4]提出了基于传播算子方法的递推预测子空间辨识方法.对于存在负载扰动的系统,文献[5]采用麦克劳林时间序列逼近带有负载扰动的输出响应,从而提出了带有偏差消除估计性能的子空间辨识方法.在文献[8]中,作者通过旋转当前的信号子空间预测将来的子空间,从而递推更新系统矩阵.文献[9]提出了一种基于可变遗忘因子的稳定递推规范变量状态空间方法.结合规范相关分析方法,文献[10]提出了带有遗忘因子的数据驱动递推子空间辨识方法.文献[11]研究了两种基于递推最小二乘的子空间辨识方法,并且分别获得带有负载扰动系统的可观测马尔可夫参数矩阵和负载扰动响应.针对线性参数变化系统(linear parameter varying,LPV),文献[12]提出了基于预测的张量回归子空间辨识方法,此方法保持了LPV马尔可夫参数矩阵的固有结构并且避免了维数灾难.

上述子空间辨识方法解决了在线递推辨识问题,它们大多研究的是系统阶数已知的离散时间系统递推更新.然而,连续时间系统的辨识与离散时间系统的辨识相比较,在参数估计、H∞控制和自校正控制等方面有更加广泛的应用.在文献[13]中,作者研究了基于拉盖尔滤波器的递推连续时间子空间辨识方法,其中采用拉盖尔滤波器解决了连续时间系统中时间导数的问题.针对连续时间系统的辨识问题,文献[14]提出了基于拉盖尔滤波器的核范数子空间辨识新方法.文献[15]采用一组泊松矩函数的线性滤波器解决了连续时间系统的导数和微分问题,结合泊松矩函数和子空间辨识方法提出了辨识新方法.文献[16]提出了基于复杂频域分析的子空间辨识方法解决了连续时间系统的辨识问题.采用随机分布理论方法可以解决系统时间导数问题,此方法具有更少的参数,并且避免了辨识过程中矩阵对数的复杂运算[6].在实际生产过程中,采用子空间辨识法建立的离线模型并不能有效准确地跟踪系统的动态变化.采用奇异值分解等线性代数工具,增加算法数值鲁棒性的同时,也相应增加了子空间辨识的在线递推困难.并且过程噪声和测量噪声的存在给连续时间系统的在线辨识带来了一定的困难.因此,当系统阶数未知时,采用递推子空间辨识方法解决连续时间系统的辨识问题是一个挑战.本文针对以上连续时间系统辨识过程中存在的具体问题,提出一种基于随机分布理论的递推子空间辨识方法.首先,采用随机分布理论重新构建连续时间系统的输入输出Hankel矩阵,解决了噪声存在带来的辨识困难以及连续时间系统的时间导数问题;然后固定数据矩阵“R”规模,递推更新数据矩阵获得系统阶数,从而降低了辨识算法的计算负担,解决了传统的奇异值分解方法带来的在线递推困难.

本文其他部分构成如下:第2节介绍问题描述;第3节介绍连续时间系统的递推子空间辨识;第4节通过仿真验证所提方法的有效性和精确性;第5节给出结论.

2 问题描述

考虑如下连续时间系统:

其中:X(t)∈Rn为状态向量,U(t)∈Rm为连续输入向量,Y(t)∈Rl为连续输出向量,A ∈Rn×n,B ∈Rn×m,C ∈Rl×n,D ∈Rl×m为系统矩阵,W(t)∈Rn,V(t)∈Rl分别为系统和观测方程的连续高斯白噪声,均值为零,协方差矩阵为

其中:E(·)是期望值算子,Q(t),S(t),R(t)是噪声强度,δ(t −τ)是δ函数.U(t)与随机噪声过程W(t),V(t)都独立.本文的目的是根据输入输出数据递推更新系统的阶数和系统矩阵A,B,C,D.在辨识连续时间系统问题之前,首先考虑以下确定性的情况:

其中输入Ud(t)是足够次可微的.考虑Yd(t)在每次不同时刻t1,t2,···,tN(不必是等距的)的直到i −1阶高阶导数,有以下输入输出矩阵关系:

通过对输入输出投影矩阵的奇异值分解可以得到矩阵Γi的列空间以及系统阶数,进而得到系统矩阵和噪声强度[7].辨识过程中,在增加算法数值鲁棒性的同时,也增加了子空间辨识的在线递推困难.辨识算法要求构建输入输出数据关系的充分条件是输入U(t)和输出Y(t)是足够次可微的,但是在讨论连续时间系统的情况下,连续高斯白噪声W(t),V(t)处处不可微,而且状态向量X(t)、输出向量Y(t)处处不可导.所以在随机噪声存在的情况下,如何设计辨识算法递推更新连续时间系统是一个难点问题.本文针对以上问题深入研究基于随机分布理论的递推子空间辨识方法(recursive subspace identification based on random distribution theory,RDT–RSID).

3 连续时间系统的递推子空间辨识

本节利用测试函数和广义分布函数理论推导出连续时间系统的输入输出矩阵方程,并采用将输入输出数据矩阵“R”规模固定的方法更新连续时间系统.具体地,首先,对连续随机分布函数进行求导,然后,得到连续时间系统的输入输出矩阵方程,最后递推更新连续时间系统.

3.1 连续随机分布函数的微分计算

定义(Ω,F,P)为概率空间,空间D是在R上的C∞-函数φ(t)所构成的空间,函数φ(t)为测试函数.一个连续随机过程{y(t,ω),−∞

如果对于任意函数φ(t)∈D,随机变量y(φ)是高斯的,即随机分布y(φ)称为高斯过程,亦为规则广义分布函数.

对于∀φ ∈D,随机过程y(t)在分布意义下的一阶二阶导数可以计算为

随机过程y(t)的k阶导数为

测试函数(φ ∈S)可以从更广泛的一类函数来选择,要求函数有足够阶的导数并且在区间(−∞,∞)上为快速递减函数[17–18].

定义{tj,j1,2,···,N}(−∞

随机过程y(t)在分布意义下的一阶、二阶导数为

随机过程y(t)的k阶导数为

3.2 输入输出矩阵方程的建立

定义{tj,j1,2,···,N}(−∞

并且X(φ)(tj),U(φ)(tj),V(φ)(tj)同样定义.

对方程(1)关于t进行微分,然后对两边结果同时乘φ(t;tj),并且在(−∞,∞)上进行积分,所以方程(1)在分布意义下的一阶导数为

代入方程(1),并结合方程(7),可以得到

结合方程(8),二阶导数以相同方式给出

结合方程(9),重复求导(i −1)次,则第(i −1)阶导数为

通过列向量来表示方程(18),可以得到输入输出矩阵方程

U(φ)(tj),W(φ)(tj)和V(φ)(tj)与Y(φ)(tj)相似定义,Γi和Hi与方程(4)相似定义,Σi定义如下:

因此,通过输出列向量表示的输入输出矩阵方程为

其中φ ∈D,状态矩阵为

输出矩阵为

并且Ui(φ),Wi(φ),Vi(φ)和Yi(φ)相似定义[19].

3.3 系统矩阵和噪声强度的递推估计

根据新的输入输出数据{ut+1,yt+1},构造数据向量yi(t+1)和ui(t+1)作为更新向量:

所以,有以下关系:

其中:ui(t −j+1)和yi(t −j+1)与yi(t+1)定义类似,u2i(t+1)和u2i(t −j+1)是辅助变量,由输入输出数据构建,

Yi(·|·),Ui(·|·)是块Hankel矩阵,U2i(·|·)是由输入输出数据组成的辅助变量矩阵.

对于式(24)左边的子矩阵在时刻t处进行LQ分解,有

对于式(24)右边的子矩阵在时刻t+1处进行类似的LQ分解,有以下关系:

根据式(32)–(36),得到

通过式(39)以及对L32(t)(t)执行奇异值分解,并依据文献[14,18]辨识过程中最小二乘法和残差分析法可得系统阶数n以及系统矩阵.通过不断添加新数据来更新矩阵L32(t)(t),从而递推更新模型(1)中系统矩阵A,B,C,D和系统阶数n.至此完成了连续时间系统(1)的递推辨识过程.

本文基于随机分布理论的递推子空间辨识(RDT–RSID)方法可归纳如下:

步骤1依据式(11)–(13),得到输入信号的一阶至k阶导数;

步骤2通过式(21),得到系统的输入输出矩阵方程;

步骤3根据式(24),构造数据矩阵等式关系;

步骤4由式(26),对新得到的输入输出向量执行LQ分解得到{Lj(t),j0,1,···,4};

步骤5由式(39),计算L31(t)执行奇异值分解L32(t)递推更新模型系统矩阵和系统阶数.

4 仿真研究

4.1 数值仿真

本节通过数值例子的仿真实验验证了提出方法RDT–RSID的有效性及比较优势.考虑如下连续时间系统[20]:

e(t)为零均值白噪声,由MATLAB加入且噪声水平为20 dB,输入信号选择伪随机数二进制序列(pseudo-random binary sequence,PRBS)信号,采样时间为Ts0.01 s,执行100次蒙特卡罗仿真.

分别采用方法文献[9]中的递推子空间辨识(recursive subspace identification,RSID)与RDT–RSID对系统(40)进行辨识,分别得到RDT–RSID与真实系统的平均伯德图以及RSID与真实系统的平均伯德图如图1–2所示.可见无论频率高低,RDT–RSID比RSID得到的模型更接近真实模型(40)的伯德图,说明所提方法RDT–RSID的辨识精度更高.

图1 RDT–RSID与真实系统的平均伯德图Fig.1 The averaged bode diagram of the RDT–RSID and true system

图2 RSID与真实系统的平均伯德图Fig.2 The averaged bode diagram of the RSID and true system

4.2 连续搅拌釜反应过程应用研究

本节采用连续搅拌釜加热器系统(continuous stirred tank heater,CSTH)的辨识数据进行仿真研究,验证提出方法的有效性和精确性.连续搅拌釜加热器是加拿大阿尔伯塔大学设计的一个小型试验装置,其基本配置如图3所示[21].搅拌釜中的冷水和热水通过加热管的蒸汽进行加热,混合的液体通过一个长管子排出.假设搅拌釜里的冷热水均匀混合,水的温度和输出流的水温相同.有3个比例积分控制器来控制冷水液位、温度和流量.冷水的液位和温度是可以被控制的,所有信号都是统一量程,为(4～20)mA.

图3 连续搅拌釜加热器的配置图Fig.3 The configuration sketch of CSTH

本实验选择液位环和温度环的设定值作为输入,搅拌釜的液位和温度作为输出.温度和液位分别受到幅值为0.5 mA和1 mA的随机二进制干扰,此干扰信号由MATLAB系统辨识工具箱中的idinput函数产生.采样时间设定为Ts1 s,采样数据长度N800.

采用方法文献[13]中的基于Laguerre滤波器的递归子空间辨识(recursive subspace identification using Laguerre filters,RSILF)和所提方法RDT–RSID进行辨识,通过RDT–RSID模型和RSILF模型的预测输出分别如图4–5所示.从图中可以看出,RDT–RSID预测数据和测量数据的轨迹具有高度的一致性.为了更方便直观的观察,给出了两种方法所对应的实际测量值和预测输出值之间的误差如图6–7所示.从图中可以看出,RDT–RSID模型输出值与测量值的误差较小,表明所提出的递推辨识方法具有较高的辨识精度.

图4 RDT–RSID模型的输出Fig.4 The outputs of the RDT–RSID

图5 RSILF模型的输出Fig.5 The outputs of the RSILF

图6 RDT–RSID模型输出值与测量值的误差Fig.6 The errors between the outputs of the RSILF and measurements

图7 RSILF模型输出值与测量值的误差Fig.7 The errors between the outputs of the RSILF and measurements

根据文献[7],辨识模型的预测误差可以通过以下方程进行计算.均方根误差(root mean squared error,RMSE)和预测误差(sum of squares of prediction error,SSPE)平方和分别如下:

基于方法RDT–RSID和基于RSILF的液位和温度测量变量的均方根误差(RMSE)和预测误差平方和(SSPE)如表1所示.

表1 RDT–RSID 和RSILF 方法建模误差的RMSE和SSPETable 1 The RMSE and SSPE of prediction error from RDT–RSID and RSILF

从表中结果可以看出,基于RDT–RSID的液位和温度测量变量的RMSE和SSPE值都比基于RSILF方法的小,这表明提出方法RDT–RSID有更好的模型辨识性能.

5 结论

本文提出了基于随机分布理论的递推子空间辨识方法.利用随机分布理论构造输入输出数据Hankel矩阵,得到分布意义下的输入输出矩阵方程,从而解决了连续时间系统的时间导数计算问题.采用将输入输出数据矩阵“R”规模固定的方法,降低了辨识算法的计算量和存储成本.通过仿真实验验证了提出方法的有效性和精确性.下一步计划在本文基础上进一步考虑输入信号和噪声的相关性对算法的影响,通过理论分析结合实验仿真说明算法的辨识效果.