1 引言

时滞在实际系统(如热工系统，化工过程等)中普遍存在，且在很多情况下未知，而时滞对系统的控制有较大影响； 另外，在系统模型已知的情况下，参数辨识是对系统进行控制的前提； 因此参数和时滞的辨识具有十分重要的实际意义^{[1, 2]}. 通常，辨识系统参数前，要预先估计出系统时滞. 但是，对于具有大时滞的多输入单输出(multiple-input single-output，MISO)系统，常用的时滞估计方法可能会有较大的近似误差^{[3, 4]}. 此外，采用常规算法辨识参数^{[5, 6, 7]}和时滞时需要大量的采样数据，增加了采样成本. 因此，对于含有较大时滞的多变量系统，为了提高辨识效率，找到一种能够在有限的采样数据下有效辨识系统参数和时滞的方法是非常必要的.

实际系统中模型的阶次通常较低，也就是说，对于含有较大时滞的多变量系统，参数化得到的高维参数向量中仅含有少量非零参数且非零参数的位置未知，这样的参数向量称为稀疏向量； 由稀疏向量表示的系统称为稀疏系统^{[8, 9, 10]}. 压缩感知(compressed sensing，CS)理论^{[11, 12, 13, 14, 15]}是研究稀疏系统的一种重要理论，CS重构理论能够根据低维的观测数据和测量矩阵重构未知的高维稀疏信号^{[16, 17]}； 在许多领域基于该理论对稀疏系统已经进行了深入的研究^{[18, 19]}，然而在控制理论的研究还不多，文^[8]针对含有时滞的确定性ARX(Auto Regressive with eXternal input)模型，在有限采样数据量下采用分块正交匹配追踪算法辨识系统参数和时滞； 文^[20]针对含有时滞的单输入单输出ARX模型，在有限的采样数据量下提出了一种l₁正则化算法辨识系统参数和时滞； 文^[21]针对噪声干扰中存在离群点的情况，采用基于压缩感知的方法对系统进行辨识.

本文针对含有未知时滞的多输入单输出有限脉冲响应系统，研究在有限的采样数据量下参数和时滞的辨识问题，提出一种梯度追踪算法，同时估计系统的参数和时滞. 文^[8]中采用分块正交匹配追踪算法估计系统的非零参数，但该算法需要计算协方差矩阵，特别是当协方差矩阵维数较大时，计算量较大^[22]. 针对这一问题，本文考虑将匹配追踪思想和梯度搜索方法相结合以避免求解协方差矩阵，从而减小计算量，提高辨识效率.

2 系统描述

考虑如下MISO-FIR系统：

其中，u_i(t)是系统第i个通道的输入，y(t)是系统的输出，d_i是第i个输入通道的时滞，v(t)是零均值方差为σ²的白噪声，B_i(z)是时间移位算子z^-1(i.e.， z^-1u(t)=u(t-1))的常系数多项式. 其中，阶次n_bi已知，参数b_ij(j=1，2，…，n_bi)以及输入通道时滞d_i未知. 不失一般性，假设当t≤0时，y(t)=0，u_i(t)=0，v(t)=0.

若不直接考虑时滞d_i，定义：

其中，l为输入数据回归长度，l的取值应大于系统最大时滞d_i和最高阶次n_bi之和，即l>max(d_i+n_bi)，则系统模型(1)可改写为 由式(2)、 式(3)可知，参数向量θ是仅有少量非零参数的高维参数向量，这样的向量是稀疏向量； 式(4)所示的系统称为稀疏系统^[8].

当有m组采样数据，即t=1，2，…，m时，定义：

可得式(4)的堆积形式：

如果m$ \gg $N，可采用最小二乘方法直接得到参数向量的估计：

然而，式(6)中参数向量θ的维数N较高，若采用最小二乘算法将耗费大量的采样时间，增加辨识成本. 本文将压缩感知重构理论应用于该系统的辨识，仅需要有限的采样数据量(甚至m<N)，即可有效地估计系统的参数和时滞.

3 梯度追踪算法

根据压缩感知重构，稀疏系统(6)的辨识问题可描述为

其中，ε>0为设定的允许误差上界. 上式中θ₀表示θ的0范数，定义为║θ║₀：=｛＃i|θ_i≠0，i=1，2，…，N｝=K，其中＃表示计个数，即0范数表示参数向量中非零元素的数目，而K称为θ的稀疏度^[10].

令φ_i表示信息矩阵Φ的第i列； θ_i表示参数向量θ的第i个参数； 式(6)中输出向量Y可改写为信息矩阵各列的线性组合加噪声向量的形式，即：

由于θ的稀疏度为K，因此上式右侧只有K个非零项. 辨识的目的就是确定上式右侧这K个非零项的位置，并根据输入输出数据估计非零项的参数，进而根据已知条件估计系统的未知时滞.

令k=1，2，…为迭代变量，λ_k为第k次迭代选出φ_i的索引，索引集Λ_k=｛λ₁，…，λ_k｝； Φ_Λk为Λ_k指示的φ_i所构成的信息矩阵； ${\hat \theta }$_k表示第k次迭代得到的参数估计； ${\hat \vartheta }$_k∈R^N是根据Λ_k的指示由${\hat \theta }$_k恢复的参数估计，${\hat \vartheta }$₀=0； 每次迭代的残差用r_k表示，初始残差r₀=Y，残差r_k由下式计算：

定义准则函数：

使ξ(i)最小的φ_iθ_i为第k次迭代选出的非零项. 极小化ξ(i)可得 将θ_i代入ξ(i)可得 由式(9)可知，若要ξ(i)最小则需残差r_k与$\frac{\phi }{{\parallel {\phi _i}\parallel }}$的内积最大，因此第k次迭代选择的非零项的位置可由下式确定： 式中，λ_k用于指示选出的φ_i的位置，并将选出的φ_i用φ_λk表示. 根据λ_k更新Λ_k： 并根据Λ_k获取子信息矩阵Φ_Λk.

定义准则函数：

应用最小二乘原理，将J(θ_k)对θ_k求导并令其为0，则有：

由式(13)可得第k次迭代的参数估计：

根据得到的${\hat \theta }$_k，由式(8)更新残差r_k. 然后按照Λ_k的指示将${\hat \theta }$_k恢复为N维的参数向量${\hat \vartheta }$_k，并同${\hat \vartheta }$_k-1比较，若║${\hat \vartheta }$_k-${\hat \vartheta }$_k-1║>ε₀(式中ε₀>0为设定的允许误差)，则令k=k+1继续迭代； 否则停止迭代，得到参数估计_k.

值得注意的是，用${\hat \theta }$_k代替式(13)中的θ_k，并将式(8)代入式(13)，则有：

上式说明残差r_k与Φ_Λk中所有的列正交； 因此，式(8)、 式(10)～式(14)的算法称为正交匹配追踪(orthogonal matching pursuit，OMP)算法. 然而，OMP算法每次迭代需要计算协方差矩阵S_k：=[Φ^T_ΛkΦ_Λk]^-1，随着k的增加，计算量将显著增加. 针对这一问题，下面采用最速下降法求解参数估计.

采用负梯度搜索原理极小化J(θ_k)得

式中，μ_k是迭代步长. 然而，上式并未用到第k次迭代选择的信息向量，且${\hat \theta }$_k-1是k-1维的并不能更新第k次迭代的参数向量； 为了解决这一问题，在更新参数估计值前，根据Λ_k在${\hat \theta }$_k-1中增加一个零元素构造一个k维的参数向量${\hat \theta }$_Λk，并用${\hat \theta }$_Λk代替式(15)中的${\hat \theta }$_k-1，用Φ_Λk代替Φ_Λk-1，则有 值得说明的是，${\hat \theta }$_Λk中与第k次选择的φ_i相对应的参数为0，即有Φ_ΛkΛk=Φ_Λk-1k-1.

为使${\hat \theta }$_k收敛，式(14)中[I-μ_kΦ^T_ΛkΦ_Λk]的所有特征值都应在单位圆内，文^[23]给出了迭代步长μ_k的范围，一种可行的保守方法是迭代步长μ_k满足：

式中，${\bar \lambda }$_max(X)表示对称矩阵X的最大特征值. 下面求解最佳迭代步长，将式(16)代入式(12)得 极小化g(μ_k)并令其导数为0，则有： 可得最佳迭代步长μ_k： 式中，$\nabla $_k：=Φ^T_Λk(Y${\bar \varphi }$_Λk${\hat \theta }$_Λk). 残差r_k可由式(8)计算，但式(8)需要进行矩阵和向量的乘法运算，计算量较大. 这里将式(16)代入式(8)，有： 上式为r_k的递推表达式，由于Φ_Λk$\nabla $_k和r_k-1已经得到，因此，式(16)同式(8)相比可减少计算量. 式(10)、 式(16)～式(19)构成了含有未知时滞MISO-FIR模型的梯度追踪(gradient pursuit，GP)算法，具体实施步骤如下：

(1) 收集有限的输入输出数据｛u_i(t)，y(t)： i=1，2，…，r； t=1，2，…，m｝，设定参数估计精度ε₀.

(2) 由式(5)构造堆积输出向量Y，信息矩阵Φ.

(3) 令k=k+1，设置初值，令r₀=Y，Λ₀=∅，${\hat \vartheta }$₀=0.

(4) 根据式(10)确定索引λ_k，并根据式(11)更新集合Λ_k，由Λ_k获取子信息矩阵Φ_Λk.

(5) 根据Λ_k由${\hat \vartheta }$_k构造${\hat \theta }$_Λk.

(6) 由式(18)计算第k次的迭代步长μ_k，由式(16)计算第k次迭代的参数估计${\hat \theta }$_k； 根据Λ_k将${\hat \theta }$_k恢复为N维的${\hat \vartheta }$_k.

(7) 比较${\hat \vartheta }$_k和${\hat \vartheta }$_k-1，如果║${\hat \vartheta }$_k-${\hat \vartheta }$_k-1║>ε₀由式(16)更新残差r_k，令k=k+1回到步骤(4)继续迭代； 否则迭代停止，获得参数估计向量${\hat \vartheta }$_k.

根据得到的参数估计值${\hat \vartheta }$∈R^N，结合输入通道的数据回归长度l以及阶次n_bi，则可估计出系统每个输入通道的时滞. 由式(2)、 式(3)可知，中共有r+1个零元素块，用n_i(i=1，2，…，r+1)表示零元素块中零元的个数，则输入通道的时滞估计为

4 性能分析

收敛性分析是辨识算法研究的一项重要内容. 本小节给出GP算法的收敛性能及其计算量分析.

定理1  对于GP算法(10)，式(16)～式(19)，在每次迭代中总存在依赖于子信息矩阵Φ_Λk的常数0<c_k<1，使得下式成立：

证明  式(19)等号两边取范数的平方可得^[24]

将式(18)代入上式得

由于$\nabla $_k=Φ^T_Λkr_k-1，根据矩阵范数的相容性，上式可写为

式中，║Φ^T_Λkr_k-1║_∞表示矩阵Φ^T_Λkr_k-1的∞范数. 由文^[25]中的定理9、10可知，对于任意的残差r_k，存在常数ω>0使得║Φ^Tr_k║_∞²>ω║r_k║²成立. 由Λ_k的选择过程可知： ║Φ^Tr_k║_∞²=║Φ^T_Λkr_k║_∞². 因此，对式(21)存在常数ω使得 令c_k=(1-ω/║Φ_Λk║²)，则定理1得证.

与OMP算法相比，GP算法每次迭代不需要计算协方差矩阵S_k=(Φ^T_ΛkΦ_Λk)^-1，因而有较小的计算量. OMP算法和GP算法的计算量如表 1所示.

5 仿真例子

考虑如下MISO-FIR系统：

系统各输入通道的时滞分别为d₁=20，d₂=3，d₃=32，d₄=47，d₅=10. l=50，系统的参数向量为

式中，0_i表示有i个零元素的行向量，θ的维数N=lr=250，稀疏度为K=$\sum {_{i = 1}^5} $n_bi=10.

仿真时，输入u_i(t)采用随机不相关的零均值单位方差信号，噪声v(t)采用均值为零方差为σ²的白噪声序列. 定义参数估计误差δ：=║θ-${\hat \vartheta }$_k║/║θ║. 取采样数据长度m=200，噪声方差为σ²=0.10²，分别采用OMP算法和GP算法辨识系统，参数估计误差δ随迭代次数k的变化如表 2所示； GP算法的参数估计误差δ随迭代次数k的变化如图 1所示.

梯度追踪算法得到的稀疏参数向量如表 3所示，最终辨识得到的参数估计为

表 3中第一行是辨识得到的稀疏参数向量中非零参数的位置，中零元素块的数目分别为n₁=20，n₂=31，n₃=77，n₄=63，n₅=11，n₆=38. 结合式(20)可得系统的时滞估计为

定义节约的计算量为η：=$\frac{{{N_1} - {N_2}}}{{{N_1}}}$×100％. 当k=15时，OMP算法和GP算法的计算量以及η分别为

综合表 1、 表 2、 图 1，式(24)和式(25)可知：

(1) GP算法和OMP算法的参数估计误差随迭代次数k的增加都逐渐减少，最终趋于稳定；

(2) 该稀疏系统在有限采样数据量下，GP算法的辨识精度接近OMP算法，见表 2；

(3) GP算法能够在有效估计系统参数的同时还能有效地估计系统的未知时滞$\hat{d}$_i；

(4) GP算法能够有效地减少计算量，见式(25).

6 结论

针对含有未知时滞的MISO-FIR系统的参数与时滞辨识问题，本文提出了一种梯度追踪算法辨识系统参数和时滞，对于大时滞的系统，该算法可以在有限的采样数据下同时估计出系统的参数与时滞. 性能分析表明该算法收敛且具有较小的计算量. 数值仿真验证了算法的有效性.