0 引言

许多工业过程, 如蒸馏塔、热交换器和反应堆等, 都可以用多变量时滞系统建模.时滞的存在使得系统的输出不能及时响应, 且稳定性变差.设计系统控制器时对时滞参数的精度要求较高, 因此, 时滞的估计和系统参数估计都是系统辨识的重要内容.近年来, 时滞与参数的联合估计方法不断涌现.文献[1]利用一种自适应滤波方法辨识系统参数和时滞; 文献[2]根据系统的阶跃响应, 使用直接辨识方法实现系统参数与时滞联合估计; 文献[3]利用递阶辨识原理, 分别采用牛顿方法和梯度方法对时滞参数和系统参数进行递阶交互估计.文献[1-3]辨识的是单输入单输出系统, 其方法推广至多变量系统时, 由于系统参数的增加, 且每个通道的时滞不一样, 算法的复杂性和计算量都将增大.因此, 多变量系统的时滞与参数的联合估计是一个具有挑战性的课题.

多变量系统可以分解为多个多输入单输出系统.对于具有未知时滞的多输入单输出系统, 若采用过参数化方法, 则辨识模型中参数向量的维数将扩大且包含许多零元, 即参数向量是稀疏的^[4-5].辨识高维模型时, 常规的辨识方法需要大量的观测数据, 使得辨识的成本增加.考虑到向量的稀疏特性, 可利用贪婪算法进行稀疏向量的估计^[6].贪婪算法包括正交匹配追踪(OMP)算法^[7]、梯度追踪算法^[8]和压缩采样匹配追踪(CoSaMP)算法^[9]等, 其特点是重构速度快且易于实现.近年来, 贪婪算法用于多输入系统的参数与时滞辨识已有一些成果.文献[10-11]对OMP算法进行改进, 分别辨识含有未知时滞的多输入单输出有限脉冲响应(MISO-FIR)系统和多输入单输出受控自回归系统; 文献[12]针对时滞MISO-FIR系统, 结合梯度搜索和匹配追踪的思想, 提出一种梯度追踪算法以获得系统参数和时滞联合估计; 文献[13]利用辅助变量方法改进CoSaMP算法, 结合模型等价原理实现一类闭环系统辨识.这类改进的贪婪算法仅需少量的采样数据便能达到辨识效果, 因而减少了辨识成本.

输出误差模型, 即系统的理想脉冲传递函数加上观测噪声, 是常用的系统模型表达形式.其辨识的难点在于, 系统的输出相对于系统特征多项式的系数是非线性的.本文基于辅助模型思想^[14-16]和OMP算法, 将最小二乘迭代方法加以改进, 研究多输入单输出输出误差(MISO-OE)系统的时滞与参数联合估计.

1 问题描述

考虑图 1所示的MISO-OE系统.

(1)

其中: u_i(t)是系统第i个通道的输入; y(t)是系统的输出; v(t)是均值为零、方差为σ²的白噪声; d_i是第i个输入通道的时滞; A_i(z)、B_i(z)是单位后移算子z^-1[i.e., u(t)z^-1=u(t-1)]的常系数时不变多项式, 即

阶次n_ai、n_bi已知, 参数a_ij、b_ij以及输入通道时滞d_i未知.不失一般性, 假设当t≤0时, y(t)=0, u_i(t)=0, v(t)=0.

辨识目标是根据测量得到的系统的输入输出数据估计系统的参数和时滞.对于系统(1), 若直接利用最小二乘法进行辨识, 可将式(1)两边进行通分(两边最多同乘以公分母A₁(z)A₂(z)... A_r(z)), 这样会大大增加系统的参数数量, 辨识计算量也将增大.为了提高辨识效率, 可采用辅助模型辨识思想^[14], 即定义中间变量

(2)

由于输入通道时滞d_i未知, 采用过参数化方法, 即令l为输入数据回归长度, l的取值应不小于各子系统时滞d_i与阶次n_bi之和的最大值, 即l≥{max(d_i+n_bi)^[11].定义子系统信息向量φ_i(t)和参数向量θ_i如下:

(3)

则式(2)可改写为

(4)

定义系统信息向量φ(t)和参数向量θ如下:

(5)

则系统(1)可写为

(6)

当有m组采样数据时, 定义堆积向量和矩阵为

(7)

则式(6)可改写为以下堆积形式:

(8)

由式(3)和(5)可知, 参数向量θ中有大量的零参数, 因此是稀疏向量, 系统(8)为稀疏系统且稀疏度K:=.由于稀疏模型(8)维数较高, 辨识计算量和难度相应增大.

2 辨识算法 2.1 辅助模型最小二乘迭代算法

假设系统模型(1)中各输入通道的时滞d_i已知, 则可采用基于辅助模型的最小二乘迭代(AM-LSI)算法估计系统参数.基于辅助模型辨识思想, 在第k次迭代时, 未知项x_i(t)用前一次迭代的辅助模型输出估计值代替, 可得AM-LSI算法^[17]如下:

若时滞未知, 则该算法在过参数化条件下, 需要花费大量的时间和精力才能获得足够的测量数据.针对稀疏模型(8), 为了提高辨识效率, 可考虑采用基于OMP算法的辨识方法进行参数和时滞估计^[10-11].然而, 由于信息矩阵中含有未知的中间变量, 利用OMP算法无法直接得到系统参数的有效估计^[18].本文基于辅助模型辨识思想, 结合OMP算法, 提出一种基于辅助模型的正交匹配追踪迭代算法.

2.2 辅助模型正交匹配追踪迭代算法

假设系统模型(1)中各输入通道的时滞d_i未知.根据压缩感知重构理论^[6], 系统(8)的辨识问题可描述为

(9)

其中: 表示θ的零范数, 为设定的允许误差.

式(9)可以采用OMP算法求解.然而, 由于中间变量x_i(t)未知, 考虑结合AM-LSI算法中的辅助模型辨识思想, 将OMP算法嵌套于迭代算法, 即通过内部OMP算法迭代得到参数估计, 并刷新中间变量x_i(t).

令k=1, 2, ...为外部迭代变量, 定义

(10)

(11)

(12)

在第k次迭代时, 用代替式(8)中的Φ, 令表示信息矩阵估计的第j列, θ_j表示参数向量θ的第j个参数, 则式(8)可改写为以下形式:

令q=1, 2, ...为内部迭代变量, λ_{k, q}为第q次迭代选出的索引, 索引集Λ_{k, q}是前q次迭代选出的索引的集合, 即Λ_{k, q}={λ_{k, 1}, λ_{k, 2}, ..., λ_{k, q}}, 为λ_{k, q}指示的所构成的子信息矩阵, 表示第q次迭代得到的参数估计, r_{k, q}表示第q次迭代的残差, 初始残差设为r_{k, 0}=Y, r_{k, q}可由下式计算:

(13)

定义准则函数

使ε(j)最小的θ_j即为第q次迭代选出的非零项.将ε(j)对θ_j求导并令其为零, 可得

将θ_j代入ε(j)中, 可得

由上式可知, 若要ε(j)最小, 则残差r_{k, q-1}与的内积需最大.因此, 第q次迭代选择的索引λ_{k, q}由下式确定:

(14)

将λ_{k, q}选出的用表示.由λ_{k, q}和分别更新索引集Λ_{k, q}和子信息矩阵, 即

(15)

(16)

定义准则函数

将J(θ_{k, q})对θ_{k, q}求导并令其为零, 可得

(17)

从而得到第q次迭代的参数估计

(18)

用代替式(17)中的θ_{k, q}, 将式(13)代入(17), 有

上式说明第q次迭代的残差r_{k, q}与中所有的列均正交, 因此在第q+1次迭代时, 只需计算r_{k, q}与中除去的各列的内积^[11].于是式(14)可改写为

(19)

当内部迭代次数q=K时, 按照Λ_{k, q}的指示将恢复为N维的参数向量, 即为外部第k次迭代的参数估计.令为第i个通道第k次迭代的参数估计, 由于, 如果用和分别代替式(4)中的φ_i(t)和θ_i, 则可以用下式刷新每个输入通道的辅助变量输出值

(20)

若则由式(20)刷新, 令k:=k+1, 继续外部迭代; 否则, 算法终止, 得到参数估计

式(7)、(10) ~ (13)和(15) ~ (20)构成辨识MISO-OE系统的辅助模型正交匹配追踪迭代(AM-OMPI)算法, 具体实施步骤如下.

Step 1:确定l, 收集有限的输入输出数据{u_i(t), y(t):i=1, 2, ..., r; t=1, 2, ..., m, 设定参数估计精度ε₀.

Step 2:由式(7)构造堆积输出向量Y.

Step 3:初始化外部迭代, 令k=1, 随机数.

Step 4:由式(10)~(12)分别构造和

Step 5:初始化内部迭代, 令q=1, r_{k, 0}=Y, Λ_{k, 0}=∅,

Step 6:由式(19)确定索引λ_{k, q}, 并由式(15)和(16)更新索引集Λ_{k, q}和子信息矩阵

Step 7:由式(18)计算第q次迭代的参数估计向量

Step 8:由式(13)更新第q次迭代的残差r_{k, q}.

Step 9:若q < K, 则令q=q+1, 回到Step 6继续内部迭代; 否则内部迭代停止, 按照Λ_{k, q}的指示将恢复为N维的

Step 10:比较和:若, 则由式(20)更新计算辅助模型输出x_ik(t), 令k=k+1, 回到Step 4继续外部迭代; 否则迭代停止, 获得参数估计向量

根据得到的参数向量估计中各零元素块的位置以及零元素的数目, 可以估计出系统未知时滞.由式(3)和(5)可知, 中共有2r个零元素块, 令零元素块中零元的个数为n_j(j=1, 2, ..., 2r), 则输入通道的时滞估计为

(21)

3 计算量分析

AM-OMPI算法完成一次外部迭代的计算量为N₁, 如表 1所示.与AM-LSI算法相比, AM-OMPI算法用内部迭代OMP算法代替一次完成最小二乘算法, 因此, AM-LSI算法完成一次迭代的计算量为由于维数N很大, 显然N₁≪N₂.

4 仿真实验

考虑如下时滞MISO-OE系统:

各输入通道的时滞分别为d₁=20, d₂=10, d₃=30.取数据回归长度l=50, 则有

其中: 表示由i个零元素构成的零元素块, 稀疏度

仿真时, 输入{u_i(t)}采用随机不相关的零均值单位方差信号, 噪声{v(t)}采用零均值方差为σ²的白噪声序列, 定义参数估计误差分别用AM-OMPI算法和AM-LSI算法辨识系统参数.当噪声方差σ²=0.10²时, 不同采样数据长度下, 两种算法的参数估计误差如表 2所示.由表 2可以看出, AM-OMPI算法的参数估计误差远小于AM-LSI算法.

当采样数据长度m=140, 噪声方差σ²=0.10²时, 两种算法的参数估计误差δ随迭代次数k的变化如图 2所示.由图 2可以看出, 采样数据有限(m < N)时, AM-LSI算法辨识精度较低, 而AM-OMPI算法只需少量的迭代便可获得良好的参数估计.另外, 根据表 1, 当m=140时, AM-OMPI算法和AM-LSI算法每次迭代的计算量分别为N₁=238 001和N₂=14 494 269, 可见, AM-OMPI算法具有计算量小的优势.

当采样数据长度m=230, 噪声方差分别为σ²=0.10²、σ²=0.50²和σ²=0.80²时, AM-OMPI算法的参数估计误差δ随迭代次数k的变化如图 3所示.经过15次迭代, 参数估计误差分别为δ=0.618 4 %, δ=3.089 3 %和δ=8.465 3 %, 可见, AM-OMPI算法对噪声敏感, 当噪声方差较小时可以获得较高的参数估计精度.

当采样数据长度m=90, 噪声方差σ²=0.10², 迭代次数k=15时, AM-OMPI算法的参数估计误差δ=0.388 8%, 辨识得到的参数估计为

(22)

由式(22)可知, 中有6个零元素块, 各零元素块中零元素的数目分别为n₁=20, n₂=28, n₃=10, n₄=38, n₅=30, n₆=18.由式(21)可得系统每个输入通道的时滞估计分别为

仿真结果表明, 针对具有时滞的MISO-OE模型, 当噪声方差较小时, AM-OMPI算法在有限的测量数据(m<N)下, 经过少量迭代即可获得精度较高的时滞与参数估计, 具有成本低、速度快、计算量小等优势.

5 结论

针对含有未知时滞的MISO-OE系统, 本文结合OMP算法和辅助模型思想, 对最小二乘迭代算法进行改进, 提出了一种辅助模型正交匹配追踪迭代算法.由仿真分析可以看出, 该算法对采样数据的需求量不大, 因而可以节约辨识成本, 且在噪声水平不高的情况下可以获得精度较高的时滞估计与参数估计.