子空间辨识包括N4SID、MOESP和CVA等经典子空间辨识方法，主要优点是利用鲁棒计算工具，如：QR分解和奇异值分解(singular value decomposition, SVD），可直接获得被辨识对象的状态空间方程^[1]，因此被广泛应用于工业过程的辨识建模^[2-3]以及控制监测^[4-6]. 在工业应用中，在线子空间辨识技术具有巨大的潜力^[7-10].

在线子空间辨识的主要难点在于如何递归更新QR和SVD. 早期为了减少计算量，通常采用近似方法来实现QR和SVD的更新问题. Verhaegen等^[11]采用SVD的一阶更新策略来实现子空间辨识的递归，但是该方法只有满足严格的噪声条件才能获得无偏估计. 近年来，子空间跟踪技术如近似投影子空间跟踪(projection approximation subspace tracking，PAST)和基于传播的方法(propagator-based method，PM)被引入在线子空间辨识方法. Lovera等^[12]将PAST方法用于递归更新SVD，实现在线子空间辨识. Mercère^[13]则将基于传播的方法用于SVD的递归实现. 杨华等^[14]构造具有遗忘因子的输入输出数据矩阵，将梯度型算法引入基于遗忘因子的状态子空间跟踪中，实现对广义能观矩阵的估计.除此之外，有学者采用直接更新QR和SVD，设计在线子空间辨识方法. Kameyama等^[15]采用bona-fide算法，直接更新QR分解中的分量，从而实现QR分解的递归更新. 谢磊等^[16]通过将Updating和Downdating操作有效结合，利用快速滑窗QR分解算法，减少了不必要的重复步骤. Alenany等^[17]引入递归约束最小二乘法，直接构建在线子空间辨识方法. 随着现代计算能力越来越强，利用Hessenberg QR等算法直接对QR和SVD更新，具有更大的工业应用潜力.

基于主成分分析的子空间辨识方法(SIMPCA)，由于采用主成分分析作为主要的子空间分解手段，能够很好地处理变量含误差(Error-In-Varibles, EIV)过程的建模和监测问题^[18-21]. 典型的基于主成分分析的子空间辨识方法一般都利用辅助变量来消除噪声项，但在主成分分析过程中会出现较大的偏差，导致系统参数得不到较为精确的估计. 为此，Wu等^[22]提出了SIMPCA-E算法，该方法先估计并消除噪声项，然后采用主成分分析，提取系统相关子空间矩阵.

本文在SIMPCA-E算法的基础上，利用Hessenberg QR算法对QR分解的R分量进行递归更新，并采用特征值分解代替SVD分解，减少辨识过程的计算量；采用遗忘因子机制，实现对时变过程的辨识.

1 问题描述

考虑如下EIV线性离散系统：

${{x}}(k \,+\, 1) \,=\, {{Ax}}(k) \,+\, {{Bu}}_0(k) \,+\, {{w}}(k),$

(1)

${{y}}(k) \,=\, {{Cx}}(k) \,+\, {{Du}}_0(k) \,+\, {{v}}(k).\;\;\;\;\;\;\;\;$

(2)

其中， ${{x}}(k) \in {{\bf R} ^n}$ 为系统的状态变量， ${{{u}}_0}(k) \in {{\bf R} ^m}$ 和 ${{y}}(k) \in {{\bf R} ^l}$ 分别为系统的输入和输出， ${{w}}(k) \in {{\bf R} ^n}$ 和 ${{v}}(k) \in {{\bf R} ^l}$ 分别为过程噪声和输出测量噪声，系统矩阵A、B、C和D为适当的维数. 此外，输入测量中还包含了白噪声：

${{u}}(k){\rm{ \,=\, }}{{{u}}_0}(k){\rm{ \,+\, }}{{\mu}} (k).$

为了辨识的统计一致性，一般考虑如下假设：

假设1：系统渐进稳定，如A的特征值在单位圆内.

假设2：(A，C)是可观的.

假设3：噪声 ${{w}}(k)$ 、 ${{v}}(k)$ 和 ${{\mu}} (k)$ 为零均值白噪声，同时与过去时刻的 ${{{u}}_0}(k)$ 统计独立.

假设4：所有噪声序列不相关.

本文旨在通过在线采集输入输出对 $\left\{ {{{u}}(k)} \right\}$ 和 $\left\{ {{{y}}(k)} \right\}$ ，利用遗忘因子机制，辨识 ${{A}}(k)$ 、 ${{B}}(k)$ 、 ${{C}}(k)$ 和 ${{D}}(k)$ . 为了简化公式表达，其中线性时变系统中的 ${ A}(k)$ 、 ${ B}(k)$ 、 ${ C}(k)$ 和 ${ D}(k)$ 还采用A、B、C和D表示.

2 基于PCA和噪声估计的在线子空间辨识 2.1 离线SIMPCA-E

过去和未来时刻的栈向量和Hankel矩阵定义如下：

${{{y}}_p}(k) \,=\, \left[ {\begin{array}{*{20}{c}} {{y}}(k \,-\, p) \\ {{y}}(k \,-\, p \,+\, 1) \\ \vdots \\{{y}}(k \,-\,) \\ \end{array}} \right] ,\;{{{y}}_f}(k) \,=\, \left[ {\begin{array}{*{20}{c}} {{y}}(k) \\ {{y}}(k \,+\, 1) \\ \vdots \\ {{y}}(k \,+\, f \,-\, 1) \\ \end{array}} \right] ,$

${{{Y}}_{p,N}} \,=\, \left[ {{{{y}}_p}(k)\;\;\;{{{y}}_p}(k \,+\, 1)\;\;\;\; \cdots \;\;\;\,{{{y}}_p}(k \,+\, N \,-\, p \,-\, f \,-\, 1)} \right],$

${{{Y}}_{f,N}} \,=\, \left[ {{{{y}}_f}(k)\;\;{{{y}}_f}(k \,+\, 1)\;\;\; \cdots \;\;\,{{{y}}_f}(k \,+\, N \,-\, f \,-\, 1)} \right].$

式中：p和f为用户自定义的过去和未来时刻，一般都大于系统阶次n. 同理定义 ${{{U}}_{p,N}}$ 、 ${{{U}}_{f,N}}$ 、 ${{{V}}_{f,N}}$ 、 ${{{W}}_{f,N}}$ 和 ${{{M}}_{f,N}}$ .

通过迭代式(1)和(2)，构建输入输出矩阵方程：

$\begin{split}{{{Y}}_{f,N}} \,=\, &{{{\varGamma}} _f}{{{X}}_{f,N}} \,+\, {{{H}}_f}{{{U}}_{f,N}} \,-\, {{{H}}_f}{{{M}}_{f,N}} \,+\, \\ &{{{G}}_f}{{{W}}_{f,N}} \,+\, {{{V}}_{f,N}}.\end{split}$

(3)

式中： ${{{X}}_{f,N}}\, = [{{x}}\;(k) \,{{x}}(k + 1)\;\cdots\;{{x}}(k + N - 1)]$ ， ${{{\varGamma}} _f}$ 为阶次为n的扩展观测矩阵：

${{{\varGamma}} _f} \,=\, {\left[ {{{{C}}^{\rm T}}\;\;{{({{CA}})}^{\rm T}}\;\; \cdots \;\;{{({{CA}}^{f \,-\, 1})}^{\rm T}}} \right]^{\rm T}}.$

定义如下三角矩阵 ${{{H}}_f}$ 和 ${{{G}}_f}$ :

${{{H}}_f} \,=\, \left[ {\begin{array}{*{20}{c}} {{D}}& {\bf 0}& \cdots & {\bf 0} \\ {{CB}}& {{D}}& \cdots & {\bf 0} \\ \vdots & \vdots & {\rm{}} & \vdots \\{{CA}}^{f - 2}{{B}}&{{CA}}^{f - 3}{{B}} &\cdots & {{D}} \\ \end{array}} \right],$

${{{G}}_f} \,=\, \left[ {\begin{array}{*{20}{c}} {\bf 0}& {\bf 0}& \cdots & {\bf 0} \\ {{C}}& {\bf 0}& \cdots & {\bf 0} \\ \vdots & \vdots & {\rm{}} & \vdots \\ {{CA}}^{f - 2}&{{CA}}^{f - 3}& \cdots &{\bf 0} \\ \end{array}} \right].$

通过选择较大的p，可得

${{{X}}_{f,N}} \,\approx\, {{{L}}_z}{{{Z}}_{p,N}}.$

(4)

式中： ${{{L}}_z}$ 如文献[1]中定义， ${{{Z}}_{p,N}} = {\left[ {\begin{array}{*{20}{c}} {{{Y}}_{p,N}^{\rm T}}&{{{U}}_{p,N}^{\rm T}} \end{array}} \right]^{\rm T}}$ .

结合式(3)和(4)，可得

${{{Y}}_{f,N}} \,=\, {{{\varGamma}} _f}{{{L}}_z}{{{Z}}_{p,N}} \,+\, {{{H}}_f}{{{U}}_{f,N}} \,+\, {{{E}}_{f,N}}.$

(5)

式中： ${{{E}}_{f,N}}{\rm{ = }}{{{G}}_f}{{{W}}_{f,N}} + {{{V}}_{f,N}} - {{{H}}_f}{{{M}}_{f,N}}$ . 噪声项 ${{{E}}_{f,N}}$ 可通过式(5)进行估计，进行如下QR分解：

$\begin{aligned}\left[ {\begin{array}{*{20}{c}} {{{Z}}_{p,N}} \\ {{{U}}_{f,N}} \\ {{{Y}}_{f,N}} \\ \end{array}} \right] \,=\, \left[ {\begin{array}{*{20}{l}} {{{R}}_{11,N}}&{}&{} \\ {{{R}}_{21,N}}&{{{R}}_{22,N}} &{} \\ {{{R}}_{31,N}}&{{{R}}_{32,N}}&{{{R}}_{33,N}} \\ \end{array}} \right]\left[ {\begin{array}{*{20}{l}} {{{Q}}_{1,N}} \\ {{{Q}}_{2,N}} \\ {{{Q}}_{3,N}} \\ \end{array}} \right].\end{aligned}$

(6)

定义：

$\begin{aligned}{{{R}}_{k,N}} \,=\, \left[ {\begin{array}{*{20}{c}} {{{R}}_{11,N}}&{}&{} \\ {{{R}}_{21,N}}&{{{R}}_{22,N}} &{} \\ {{{R}}_{31,N}}&{{{R}}_{32,N}}&{{{R}}_{33,N}} \\ \end{array}} \right] .\end{aligned}$

通过QR分解的结果,可得噪声项 ${{{E}}_{f,N}}$ 的估计如下：

${\hat {{E}}_{f,N}} \,=\, {{{R}}_{33,N}}{{{Q}}_{3,N}}.$

(7)

将 ${\hat {{E}}_{f,N}}$ 代入式(5)并移至等式左边：

${{{Y}}_{f,N}} \,-\, {\hat {{E}}_{f,N}} \,=\, {{{\varGamma}} _f}{{{X}}_{f,N}} \,+\, {{{H}}_f}{{{U}}_{f,N}}.$

(8)

定义：

$\begin{aligned}&{\tilde {{Y}}_{f,N}} \,=\, {{{Y}}_{f,N}} - {\hat {{E}}_{f,N}},\\&\left[ {{{I}} \,-\, {{{H}}_f}} \right]\left[ {\begin{array}{*{20}{c}} {{\tilde Y}_{f,N}} \\ {{{U}}_{f,N}} \\ \end{array}} \right] \,=\, {{{\varGamma}} _{f,N}}{{{X}}_{f,N}}.\end{aligned}$

(9)

式(9)两边乘以 ${{{\varGamma}} _f}$ 的补 ${({{\varGamma}} _{_f}^ \bot )^{\rm T}}$ ：

${({{\varGamma}} _f^ \bot )^{\rm T}}\left[ {{{I}} \,-\, {{{H}}_f}} \right]\left[ {\begin{array}{*{20}{c}} {{\tilde Y}_{f,N}} \\ {{{U}}_{f,N}} \\ \end{array}} \right] = 0.$

(10)

进行如下主成分分析：

$\left[ {\begin{array}{*{20}{c}} {{\tilde Y}_{f,N}} \\ {{{U}}_{f,N}} \\ \end{array}} \right] \,=\, {{PT}}^{\rm T} \,+\, \overline {{P}}\,{\overline {{T}}^{\rm T}}.$

(11)

式中： $\overline {{P}}$ 和 $\overline {{T}}$ 分别为残差空间的负荷矩阵和得分矩阵. 利用所估计的 $\overline {{P}}$ ，可得 ${{{\varGamma}} _f}$ 和 ${{{H}}_f}$ 的估计：

$\left[ {\begin{array}{*{20}{c}} {{\varGamma}} _{_f}^ \bot \\ - {({{{H}}_f})^T}{{\varGamma}} _{_f}^ \bot \\ \end{array}} \right] \,=\, \overline {{P}} {{F}}{\rm{ \,=\, }}\left[ {\begin{array}{*{20}{l}} \overline {{{{P}}_\Gamma }} \\ \overline {{{{P}}_U}} \\ \end{array}} \right] .$

(12)

式中：F为任意非奇异矩阵.一般选取F为单位矩阵，系统矩阵可与参考文献[23]类似，从 $\overline {{{{P}}_\Gamma }} $ 和 $\overline {{{{P}}_U}} $ 提取.本文仿真主要源于文献[22-23],对于系统阶次的确定，可在初始辨识时采用离线SIMPCA-E方法. 具体来讲，主要是利用式(11)进行奇异值分解，通过对奇异值大小进行排序，选取其中较大的奇异值个数作为系统的阶次，具体可参照文献[22]，其他方法包括赤池准则(Akaike information criterion, AIC)，误差重构方法(reconstruction error method)，可参考文献[24]. 这里为了省略，均假设该系统阶次已知.

2.2 在线SIMPCA-E 2.2.1 QR分解与特征值分解的集成

离线SIMPCA-E算法的关键是，如何估计噪声项 ${{{E}}_{f,N}}$ 同式(5)和残差子空间 $\overline {{P}}$ . 通过采用特征值分解代替SVD分解，可将噪声项和残差子空间的估计相结合.

在QR分解中，可得如下等式：

${{{U}}_{f,N}}{\rm{ \,=\, }}{{{R}}_{{\rm{21,}}N}}{{{Q}}_{1,N}} \,+\, {{{R}}_{{\rm{22,}}N}}{{{Q}}_{2,N}},$

(13)

${{{Y}}_{f,N}}{\rm{ \,=\, }}{{{R}}_{{\rm{31,}}N}}{{{Q}}_{1,N}} \,+\, {{{R}}_{{\rm{32,}}N}}{{{Q}}_{2,N}} \,+\, {{{R}}_{{\rm{33,}}N}}{{{Q}}_{3,N}}.$

(14)

因此，可构建如下矩阵：

$\begin{align} &\left[ {\begin{array}{*{20}{l}} {{{Y}}_{f,N}} \,-\, {{\hat {{E}}}_{f,N}} \\ {{{U}}_{f,N}} \\ \end{array}} \right] \left[ {{{({{{Y}}_{f,N}} \,-\, {{\hat {{E}}}_{f,N}})}^{\rm T}} {{U}}_{f,N}^{\rm T}} \right]\,=\, \\& \left[ {\begin{array}{*{20}{l}} {({{{Y}}_{f,N}} \,-\, {{\hat {{E}}}_{f,N}}){{({{{Y}}_{f,N}} \,-\, {{\hat {{E}}}_{f,N}})}^{\rm T}}}&{({{{Y}}_{f,N}} \,-\, {{\hat {{E}}}_{f,N}}){{U}}_{f,N}^{\rm T}} \\ {{{{U}}_{f,N}}{{({{{Y}}_{f,N}} \,-\, {{\hat {{E}}}_{f,N}})}^{\rm T}}}&{{{{U}}_{f,N}}{{U}}_{f,N}^{\rm T}} \end{array}} \right]. \\ \end{align} $

(15)

通过替代QR分解中的R矩阵，可得

$\begin{split}&({{{Y}}_{f,N}} \,-\, {\hat {{E}}_{f,N}}){({{{Y}}_{f,N}} \,-\, {\hat {{E}}_{f,N}})^{\rm T}} \,=\,\\ &{{{R}}_{31,N}}{{R}}_{31,N}^{\rm T} \,+\, {{{R}}_{32,N}}{{R}}_{32,N}^{\rm T},\end{split}$

(16)

$({{{Y}}_{f,N}} \,-\, {\hat {{E}}_{f,N}}){{{U}}_{f,N}}^{\rm T} \,=\, {{{R}}_{31,N}}{{R}}_{21,N}^{\rm T} \,+\, {{{R}}_{32,N}}{{R}}_{22,N}^{\rm T},$

(17)

${{{U}}_{f,N}}{({{{Y}}_{f,N}} \,-\, {\hat {{E}}_{f,N}})^{\rm T}} \,=\, {{{R}}_{21,N}}{{R}}_{31,N}^{\rm T} \,+\, {{{R}}_{22,N}}{{R}}_{32,N}^{\rm T},$

(18)

${{{U}}_{f,N}}{{U}}_{f,N}^{\rm T} \,=\, {{{R}}_{21,N}}{{R}}_{21,N}^{\rm T} \,+\, {{{R}}_{22,N}}{{R}}_{22,N}^{\rm T}.$

(19)

式(15)可转换为

$\begin{aligned} &\left[ \!\!\!\! {\begin{array}{*{20}{l}} {{{Y}}_{f,N}} \,-\, {{\hat {{E}}}_{f,N}} \\ {{{U}}_{f,N}} \\ \end{array}} \right]\left[ {{{({{{Y}}_{f,N}} \,-\, {{\hat {{E}}}_{f,N}})}^{\rm T}} {{U}}_{f,N}^{\rm T}} \right] \,=\, \\ &\left[ \!\!\!\!{\begin{array}{*{20}{l}} {{{{R}}_{31,N}}{{R}}_{31,N}^{\rm T} \,+\, {{{R}}_{32,N}}{{R}}_{32,N}^{\rm T}}\quad\;\;\; \!{{{{R}}_{31,N}}{{R}}_{21,N}^{\rm T} \,+\, {{{R}}_{32,N}}{{R}}_{22,N}^{\rm T}} \\ {{{({{{R}}_{31,N}}{{R}}_{21,N}^{\rm T} \,+\, {{{R}}_{32,N}}{{R}}_{22,N}^{\rm T})}^{\rm T}}}\;\;\! {{{{R}}_{21,N}}{{R}}_{21,N}^{\rm T} \,+\, {{{R}}_{22,N}}{{R}}_{22,N}^{\rm T}}\end{array}}\!\!\!\! \right]. \\\end{aligned} $

(20)

进行如下特征值分解：

$\left[ {\begin{array}{*{20}{c}} {{\tilde Y}_{f,N}} \\ {{{U}}_{f,N}} \\ \end{array}} \right] \left[ {{{\tilde Y}_{f,N}}^{\rm T} \quad {{U}}_{f,N}^{\rm T}} \right] \,=\, [{ \varPi _{1}}\quad { \varPi _{2}}]\left[ {\begin{array}{*{20}{c}} {{{{\varLambda}} _1}}&{\bf 0} \\ {\bf 0}&{{{{\varLambda}} _2}} \end{array}} \right]\left[ {\begin{array}{*{20}{c}} {{{\tilde { \varPi}}_1}} \\ {{{\tilde { \varPi}}_2}} \end{array}} \right].$

(21)

式中： ${{ \varPi} _1}$ 和 ${{ \varPi} _2}$ 为对应的特征值向量矩阵. ${{ \varPi} _2}$ 可表示残差子空间，因此

$\left[ {\begin{array}{*{20}{l}} {{\varGamma}} _{_f}^ \bot \\ - {({{{H}}_f})^T}{{\varGamma}} _{_f}^ \bot \\ \end{array}} \right] \,=\, { \varPi _2}{{F}}{\rm{ \,=\, }}\left[ {\begin{array}{*{20}{l}} { \varPi} _{2\varGamma } \\ { \varPi} _{2U} \\ \end{array}} \right] .$

(22)

2.2.2 在线残差子空间更新

从2.2.1节可知，QR分解的更新是在线实现SIMPCA-E算法的关键是 ${{{R}}_{31,N}}\;{{R}}_{31,N}^{\rm T}{\text{、}}$ ${{{R}}_{32,N}}\;{{R}}_{32,N}^{\rm T}{\text{、}}$ ${{{R}}_{31,N}}\;{{R}}_{21,N}^{\rm T}{\text{、}}$ ${{{R}}_{22,N}}\;{{R}}_{32,N}^{\rm T,}{\text{、}}$ ${{{R}}_{21,N}}\;{{R}}_{21,N}^{\rm T}$ 和 ${{{R}}_{22,N}}\;{{R}}_{22,N}^{\rm T}$ 的递归更新.

假设已采集初始数据 $\{ {{u}}(k)\} _{k = 1}^N$ 和 $\{ {{y}}(k)\} _{k = 1}^N$ ，由式(6)和式(16)~(20)，可以得到R矩阵. 这里，采用Hessenberg QR方法对R分量矩阵进行递归更新.

当获取了新的输入输出对 $u(N + 1)$ 和 $y(N + 1)$ ,式(6)可重构为

$\begin{aligned}&\left[ \!\!\!\! \begin{array}{l}{{ Z}_{p,N}} \quad { z_p}(k \,+\, N \,-\, p \,-\, f)\\{{ U}_{f,N}} \quad {u_f}(k \,+\, N \,-\, f) \\{{ Y}_{f,N}} \quad {y_f}(k \,+\, N \,-\, f)\end{array} \!\!\!\! \right]\,=\,\\& \left[ \!\!\!\! \begin{array}{l}{{ {z}}_p}(k \,+\, N \,-\, p \,-\, f) \,\,\, \,\,\,\, {{ R}_{11,N}}\;\;\;\\{{ u}_f}(k \,+\, N \,-\, f) \quad \quad {{ R}_{21,N}}\;\;\;\; {{ R}_{22,N}}\\{{ y}_f}(k \,+\, N \,-\, f) \quad \quad {{ R}_{31,N}}\;\;\;\;\;{{ R}_{32,N\;\;\;\;\;\;\;\;\;\;}}{{ R}_{33,N}}\end{array} \!\!\!\! \right]\left[ \!\!\!\! \begin{array}{l}{\bf{0}} \quad \,\,\,\,\,\, {\bf 1} \\{{ \bar { Q}}_1} \quad {\bf{0}}\\{{ \bar { Q}}_2} \quad {\bf{0}}\end{array} \!\!\!\! \right].\end{aligned}$

(23)

式中：

${{{z}}_p}(k \,+\, N) \,=\, \left[ {\begin{array}{*{20}{c}} {{{y}}_p}(k \,+\, N \,-\, p \,-\, f) \\ {{{u}}_p}(k \,+\, N \,-\, p \,-\, f) \\ \end{array}} \right] .$

定义：

$\begin{aligned}&{{{H}}_{k,N \,+\, 1}} \,=\,\\ &\left[{\begin{array}{*{20}{l}}{{{z}}_p}(k \,+\, N \,-\, p \,-\, f){}& {{{R}}_{11,N}}&{}&{} \\ {{{u}}_f}(k \,+\, N \,-\, f) &{{{R}}_{21,N}}&{{{R}}_{22,N}}&{} \\ {{{y}}_f}(k \,+\, N \,-\, f)& {{{R}}_{31,N}}&{{{R}}_{32,N}}&{{{R}}_{33,N}} \\ \end{array}} \right] .\end{aligned}$

(24)

R矩阵可以利用式(23)通过Hessenberg QR法进行更新，即利用一系列Givens旋转 ${q_i}$ ，可得

$\begin{aligned}&{{{H}}_{k,N + 1}}\times{{{q}}_1}{{{q}}_2} \cdots {{{q}}_x} \,=\, \\&\left[ {\begin{array}{*{20}{c}} {{{R}}_{11,N \,+\, 1}}&{}&{}& {\bf 0} \\ {{{R}}_{21,N \,+\, 1}}&{{{R}}_{22,N \,+\, 1}} &{}&{\bf 0} \\ {{{R}}_{31,N \,+\, 1}}&{{{R}}_{32,N \,+\, 1}}&{{{R}}_{33,N \,+\, 1}} &{\bf 0} \\ \end{array}} \right] .\end{aligned}$

(25)

去除最后一列，设置

${{{R}}_{k,N{\rm{ \,+\, 1}}}} \,=\, \left[ {\begin{array}{*{20}{c}} {{{R}}_{11,N \,+\, 1}}&{}&{} \\ {{{R}}_{21,N \,+\, 1}}&{{{R}}_{22,N \,+\, 1}} &{} \\ {{{R}}_{31,N \,+\, 1}}&{{{R}}_{32,N \,+\, 1}}&{{{R}}_{33,N \,+\, 1}} \\ \end{array}} \right] .$

由此可得R矩阵的更新：

$\begin{split}&({{{Y}}_{f,N \,+\, 1}} - {\hat {{E}}_{f,N \,+\, 1}}){({{{Y}}_{f,N \,+\, 1}} - {\hat {{E}}_{f,N \,+\, 1}})^{\rm T}} =\\& {{{R}}_{31,N \,+\, 1}}{{R}}_{31,N \,+\, 1}^{\rm T} \,+\, {{{R}}_{32,N{\rm{ \,+\, 1}}}}{{R}}_{32,N \,+\, 1}^{\rm T},\end{split}$

(26)

$\begin{split}&({{{Y}}_{f,N \,+\, 1}} - {\hat {{E}}_{f,N \,+\, 1}}){{{U}}_{f,N \,+\, 1}}^{\rm T}= \\& {{{R}}_{31,N \,+\, 1}}{{R}}_{21,N{\rm{ \,+\, 1}}}^{\rm T} + {{{R}}_{32,N{\rm{ + 1}}}}{{R}}_{22,N \,+\, 1}^{\rm T},\end{split}$

(27)

$\begin{split}&{{{U}}_{f,N \,+\, 1}}{({{{Y}}_{f,N \,+\, 1}} \,-\, {\hat {{E}}_{f,N \,+\, 1}})^{\rm T}}\,=\, \\& {{{R}}_{21,N \,+\, 1}}{{R}}_{31,N \,+\, 1}^{\rm T} \,+\, {{{R}}_{22,N \,+\, 1}}{{R}}_{32,N{\rm{ \,+\, 1}}}^{\rm T},\end{split}$

(28)

$\begin{split}&{{{U}}_{f,N \,+\, 1}}{{U}}_{f,N \,+\, 1}^{\rm T} \,=\,\\ &{{{R}}_{21,N \,+\, 1}}{{R}}_{21,N \,+\, 1}^{\rm T} \,+\, {{{R}}_{22,N \,+\, 1}}{{R}}_{22,N \,+\, 1}^{\rm T}.\end{split}$

(29)

当R矩阵更新后,可从式(20)~(22)获取扩展观测矩阵.

2.2.3 主要结果

综上所示，本文提出的在线SIMPCA-E算法可归纳为如下步骤.

算法1(OSIMPCA-E)：

1)利用初始数据序列 $\{ {{u}}(k)\} _{k = 1}^{{N_0}}$ 和 $\{ {{y}}(k)\} _{k = 1}^{{N_0}}$ ，通过式(6)获取 ${{{R}}_{k,{N_0}}}$ 矩阵，构建 ${{{H}}_{k,{N_0}}}$ .

2)当采集到最新的输入输出 ${{u}}({N_0} + 1)$ 和 ${{y}}({N_0} + 1)$ ，构造 ${{{H}}_{k,{N_0} + 1}}$ ，通过式(25)的Givens旋转，利用Hessenberg QR算法如式(26) ~ (29)更新 ${{{R}}_{k,{N_0}{\rm{ + 1}}}}$ 矩阵.

3)如式(21)所示，应用特征值分解，估计残差子空间 ${\Pi _2}$ .

4)从 ${\bf \varPi _2}$ 中获取 ${{{\varGamma}} _f}$ 和 ${{{H}}_f}$ ，提取A、B、C、D系统矩阵.

5)置 ${{{R}}_{k,{N_0} + 1}} \to {{{R}}_{k,{N_0}}},{N_0} + 2 \to {N_0} + 1$ ，回到步骤2).

本文所提方法的计算量主要包括QR更新和特征值分解，见(21)和(25). 当f和p相同时，Hessenberg QR算法的计算量为 $3{(2l + 2m)^2}{f^2}$ ，特征值分解的计算量为 ${(l + m)^3}{f^3}$ .

2.2.4 时变系统的跟踪

为了能够跟踪时变过程，类似于文献[23]，本文将遗忘因子机制引入所提的OSIMPCA-E算法，直接代替式(24)：

$\begin{split}&{{{H}}_{k,N \,+\, 1}} =\\ &\left[ {\begin{array}{*{20}{l}} {{{z}}_p}(k \,+\, N \,-\, p \,-\, f) &\lambda {{{R}}_{11,N}}&{}&{} \\ {{{u}}_f}(k \,+\, N \,-\, f) &\lambda {{{R}}_{21,N}}&\lambda {{{R}}_{22,N}} &{} \\ {{{y}}_f}(k \,+\, N \,-\, f) &\lambda {{{R}}_{31,N}}&\lambda {{{R}}_{32,N}}&\lambda {{{R}}_{33,N}} \\ \end{array}} \right] .\end{split}$

(30)

式中： $\lambda (0 < \lambda < 1)$ 为遗忘因子.式(30)代表最新的数据信息.

遗忘因子 $\lambda $ 通常接近于1. 一般来讲，其值采用试凑法获得^[13]. 所提OSIMPCA-E算法的关键是递归更新 ${ R}$ 矩阵，因此在引入遗忘因子后，可对其收敛性进行简单分析，以 ${{{R}}_{31,N + 1}}{{R}}_{31,N + 1}^{\rm T}$ 的更新为例，其他R矩阵类似.

假设

${{H}}_{_{k,N}}^{\rm T} \,=\, \left[ \begin{array}{*{20}{c}} {h_{11,1}} &\cdots &{h_{11,N}} \\ & \ddots & \vdots \\& {} & {}{h_{nn,N}} \end{array} \right],$

则可得

$\overline {{H}}_{_{k,N \,+\, 1}}^{\rm T} \,=\, \left[ {\begin{array}{*{20}{c}} \lambda {h_{11,1}}& \cdots &\lambda {h_{11,N}} \\& \ddots & \vdots \\&{}& \lambda {h_{nn,N}} \\ {x_1} & \cdots &{x_n} \\ \end{array}} \right] .$

经过Givens转换后，可得

${{H}}_{_{k,N + 1}}^{\rm T} \,=\, \overline {{H}}_{_{k,N \,+\, 1}}^{\rm T}{{{Q}}^{\rm T}} \,=\, \left[ {\begin{array}{*{20}{c}} {x_1}& \cdots & {x_n} \\ \lambda {h_{11,1}}& \cdots&\lambda {h_{11,N}} \\&\ddots & \vdots \\&{}& \lambda {h_{nn,N}} \\ 0 &\cdots &0 \\ \end{array}} \right] $

由 ${{H}}_{_{k,N + 1}}^{\rm T}$ 推导得

${{{R}}_{31,N \,+\, 1}}{{R}}_{31,N + 1}^{\rm T} \,=\, {\lambda ^2}{{{R}}_{31,N}}{{R}}_{31,N}^{\rm T} \,+\, {{{\varDelta}} _{N \,+\, 1}}.$

这里 ${{{\varDelta}} _{N + 1}}$ 表示与最新输入输出的相关项，而与 ${{{R}}_{31,N}}{{R}}_{31,N}^{\rm T}$ 无关. 以此类推，若遗忘因子 $0 < \lambda < 1$ ，则以前时刻 ${{{R}}_{31,N0}}$ 相关项将指数级的逐渐趋向于0，从而使得当前 ${{{R}}_{31,N + 1}}$ 相关项跟踪到系统的变化.

3 仿真结果 3.1 数值仿真

为了表明所提方法的有效性，采用文献[24]的模型，系统描述如下：

$\begin{aligned}{{x}}(t + 1) =& \left[ {\begin{array}{*{20}{c}}& 0.8 \;&- 0.4 \;&0.{\rm{2}} \\ &0&0.{\rm{3}} &- 0.5 \\ &0&0& 0.5 \\ \end{array}} \right] {{x}}(t) + \\ &\left[ {\begin{array}{*{20}{c}} &0 &0 \\ &0 &- 0.6 \\& 0.5 &0 \\ \end{array}} \right] {{u}}(t) +\left[ {\begin{array}{*{20}{c}} 0 .055 \\ 0.040 \\ 0.045 \\ \end{array}} \right] {{w}}(t),\end{aligned}$

$ {{y}}(t) = \left[ {\begin{array}{*{20}{l}} 0.5 \;\;\;0.5\;\;\;\;\;0 \\ 0\;\;\;\;\;\;0\;\;\;\;\;\;\;\;\;1 \\ \end{array}} \right] {{x}}(t) + \left[ {\begin{array}{*{20}{l}} 0.025 \\ 0.030 \\ \end{array}} \right] {{v}}(t).$

为了能够评估OSIMPCA-E算法的性能，考虑系统矩阵A不变化、缓慢变化、突变等3种情况. 假设已知系统阶次，为3. 将过去时刻p和未来时刻f设置为6. 针对每种情况，共做100 蒙特卡洛实验. 每次实验采集2 000个数据. 初始数据长度为200. 噪声信号设定如下： ${{w}}(k)$ 和 ${{v}}(k)$ 为高斯白噪声，方差为1： ${\sigma _w} = 1$ 和 ${\sigma _v} = 1$ . 输入激励信号 ${{{u}}_0}(k)$ 为单位高斯白噪声. 输入噪声 ${{\mu}} (k)$ 为高斯白噪声，方差 ${\sigma _\mu } = 0.1$ .

1) 系统矩阵A不变. 首先，针对线性时不变系统测试所提方法，选择 $\lambda = 0.998$ . 图1表示了针对时不变过程，OSIMPCA-E算法所辨识的系统极点(p_1、p_2、p_3) 与真实系统极点的对比. 图中，N_sample为样本数；p₁、p₂、p₃为系统特征平均值. 从图1可知，针对时不变过程，OSIMPCA-E算法能够获得很好的系统矩阵辨识精度. 对于大部分基于子空间辨识的过程监测来讲，由 ${\varGamma _f}$ 组成的子空间估计至关重要，因为很多过程监测方法都是采用估计的残差子空间来构建Luenberger观测器. 因此，仿真实验也对辨识所得的扩展观测矩阵与真实扩展观测矩阵所形成的子空间夹角(cos ϕ)进行计算，如图2所示. 当辨识的扩展观测矩阵与真实扩展观测矩阵相一致时，其夹角为0. 由图2可知，OSIMPCA-E算法能够获得很好的扩展观测矩阵辨识精度.

2) 系统矩阵A缓慢变化

当采样1 000个数据后, 系统矩阵A如下：

$\begin{aligned}&{{A}} =\left[ {\begin{array}{*{20}{c}} 0{\rm{.8}} - 0.2 M & - 0.4 & 0.{\rm{2}} \\ 0& 0.{\rm{3}} - 0.2 M & - 0.5 \\ 0& 0& 0.5 + 0.1 M \\ \end{array}} \right] .\end{aligned}$

式中：

$\begin{aligned}&{\displaystyle{M={{{\rm \exp} \left.\left( - \frac{{t - 1\,000}}{{2\,000}}\right)\right/}}{{{{\rm e}^{ - 1}} - 1}}}}\end{aligned}$

这里选择λ= 0.920. 图3表示了针对缓慢时变过程，OSIMPCA-E算法所辨识的系统极点与真实系统极点的对比。从图3可知，虽然在精度上还有待提高，OSIMPCA-E算法针对缓慢时变过程，能够较快的跟踪系统矩阵的变化。同时，所辨识的扩展观测矩阵也能够匹配真实扩展观测矩阵的变化，如图4所示。

3) 系统矩阵A突变

当采样1 000个数据后，系统矩阵A如下：

${{A}} = \left[ {\begin{array}{*{20}{l}} 0{\rm{.9}} & - 0.4 &0.{\rm{2}} \\0&0.{\rm{1}}& - 0.5\\ 0&0&0.{\rm{7}} \\ \end{array}} \right] .$

这里选择 $\lambda = 0.9300$ . 针对突变过程，所提OSIMPCA-E算法能够快速跟踪系统的变化，包括系统矩阵以及扩展观测矩阵，如图5和6所示. 综上所述，OSIMPCA-E针对时不变系统具有较好的辨识精度，而针对时变系统包括缓慢时变和突变过程，OSIMPCA-E算法能够很好的跟踪到系统参数的变化.

3.2 工业仿真

本文将考虑一个工业Benchmark模型来表明所提算法的有效性. 该模型为污水处理模型所采用的过程为巴西圣保罗大学所建(Activated Sludge Wastewater Treatment Plant-University of Sao Paulo，ASWWTP-USP)^[25]. 该过程的辨识目标是建立一定工况下厌氧区和好氧区的氮含量. 本文采用伪随机二进制序列信号(Pseudo-random binary sequences，PRBS)作为激励输入. 共产生1 400组数据，预处理后的数据如图7所示. 采用前300组数据进行初始建模，系统阶次选为3，过去时刻和未来时刻选为6， $\lambda = 0.999\,8$ .

图8显示了在线辨识的后1 000个数据. 结果表明OSIMPCA-E方法能够获得良好的性能，可以快速跟踪过程的变化. 为了进一步比较所提方法的性能，采用相对平均误差指标 ${\eta _{{\rm MRSE}}}$ ，其定义如下:

${\eta _{{\rm{MRSE}}}} = \frac{1}{l}{\sum_{i = 1}^l} {\;\;{{\left\{ {\frac{{{\displaystyle\sum\nolimits_{j = 1}^N} {\;{{\left[ {{y_i}(j) - {{\hat y}_i}(j)} \right]}^2}} }}{{\displaystyle\sum\nolimits_{j = 1}^N {\;{y_i}^2(j)} }}} \right\}}^{1/2}}} .$

式中：y_i表示真实输出， ${\hat y_i}$ 表示预测. 如表1所示，通过后1 000个数据进行验证，与离线算法MOESP、N4SID等^[25]比较，可发现所提算法能够获得令人满意的结果.

4 结　语

本文提出了一种基于主成分分析和噪声估计的在线子空间辨识方法. 为了减少计算量，该方法采用Hessenberg QR方法更新QR，通过主成分分析实现残差子空间的在线估计，并利用遗忘因子机制实现时变系统的辨识. 仿真结果表明，所提方法能够很好地估计时不变、时变以及突变线性EIV系统. 未来研究方向将包括基于在线子空间辨识方法构建时变过程的故障检测与诊断.