丁锋; 徐玲; 刘喜梅
【期刊名称】《《青岛科技大学学报（自然科学版）》》 【年(卷),期】2019(040)005 【总页数】20页(P1-20)
【关键词】传递函数; 参数估计; 迭代辨识; 频率特性; 梯度搜索; 多新息辨识; 耦合辨识
【作 者】丁锋; 徐玲; 刘喜梅
【作者单位】江南大学 物联网工程学院 江苏 无锡 214122; 青岛科技大学 自动化与电子工程学院 山东 青岛 266061 【正文语种】中 文 【中图分类】TP273
基于系统输入输出离散观测数据的统计辨识方法得到很大发展。为了解决存在不可测变量的系统辨识问题，以及存在有色噪声干扰的系统辨识问题等而提出了辅助模型辨识方法、多新息辨识方法、递阶辨识方法、耦合辨识方法、滤波辨识方法等[1-4]。这些方法可用于研究信号模型和传递函数的参数估计[5-7]。
在线性系统的频域分析法中，频率特性是一个重要的分析工具。控制系统的频率特性可以通过频率响应实验获得，并且可以表示为多种形式，如绘制波特图(Bode plot)、奈奎斯特图(Nyquist plot)、尼克尔斯图(Nichols plot)等。频域法分析控
制系统性能的优势在于无需求解系统的微分方程，只需要通过绘制出频率特性曲线，采用图形方式，通过频域和时域之间的关系实现系统性能直观形象的分析。 在工程实践中，往往并不需要准确地计算系统响应的全部过程，而是希望避开复杂的计算，简单、直观地分析出系统结构、参数对系统性能的影响。借助频率特性分析仪等测试手段直接获得系统的频率特性，建立数学模型作为分析与设计系统的依据，这对难以用理论分析方法去建立数学模型的系统极为有利。频率特性是一种非参数模型，控制系统的频率特性可以用实验方法测定，对于系统辨识而言，对于难以用机理法建立数学模型的系统，采用频率特性观测数据建立系统的数学模型也是一个有效方法。
连载论文[7]基于系统的频率响应，即频率特性——实频特性和虚频特性观测数据，讨论了一般线性时不变系统传递函数参数的递推辨识方法，本研究就以频率响应观测数据研究传递函数参数的迭代辨识方法。相关传递函数参数辨识方法和信号模型参数估计方法等可参见文献[8-22]。 1 一般线性定常系统的频率特性 设线性定常系统对应的传递函数模型为 (1)
其中，a1，a2，…，an，b0，b1，…，bn,为待拟合模型的参数，n为模型阶次。为分析方便，这里假设传递函数分子和分母多项式阶次相同,实际上也可以讨论阶次不同的情形。
当n为奇数时，定义参数向量和信息向量: θ1=[a1,a3,a5,…,an]T∈(n+1)/2, θ2=[a2,a4,a6,…,an-1]T∈(n-1)/2, ρ1=[b1,b3,b5,…,bn]T∈(n+1)/2,
ρ2=[b2,b4,b6,…,bn-1]T∈(n-1)/2, φ(ω)=
[ω,-ω3,ω5,-ω7,…,(-1)(n-1)/2ωn]T∈(n+1)/2, (2) ψ(ω)=
[-ω2,ω4,-ω6,…,(-1)(n-1)/2ωn-1]T∈(n-1)/2。 (3)
当n为偶数时,定义参数向量和信息向量如下: θ1=[a1,a3,a5,…,an-1]T∈n/2, θ2=[a2,a4,a6,…,an]T∈n/2, ρ1=[b1,b3,b5,…,bn-1]T∈n/2, ρ2=[b2,b4,b6,…,bn]T∈n/2, φ(ω)=
[ω,-ω3,ω5,-ω7,…,(-1)n/2-1ωn-1]T∈n/2, (4) ψ(ω)=
[-ω2,ω4,-ω6,…,(-1)n/2ωn]T∈n/2。 (5)
参考文献[7]，将s=jω代入式(1)得到模型的频率特性为
Re[G(jω)]+jIm[G(jω)], (6)
其中Re[G(jω)]和Im[G(jω)]为传递函数模型的频率特性实部和虚部，即实频特性
和虚频特性。
对于实际系统，通过频率特性试验,可以收集到每一个离散频率点ωk对应的频率特性观测数据{A(ωk)ejη(ωk)，k=0,1,2,…},或实频和虚频特性数据{M(ωk),N(ωk)，k=0,1,2,…}，且
H(jωk)=A(ωk)ejη(ωk)=M(ωk)+jN(ωk), M(ωk)=A(ωk)cosη(ωk), N(ωk)=A(ωk)sinη(ωk)。
也就是说，可以根据收集的频率特性幅值A(ωk)和相位η(ωk)计算出对应的实部和虚部数据M(ωk)及N(ωk)。下面基于频率特性数据M(ωk)+jN(ωk)，构造关于系统模型参数的准则函数，利用某种优化策略，研究估计系统传递函数G(s)参数al和bl的迭代辨识方法。 2 联立实部虚部迭代辨识方法
频率响应辨识就是找到一个传递函数模型，使其实频特性Re[G(jωk)]和虚频特性Im[G(jωk)]分别与系统的实频特性观测数据M(ωk)和虚频特性观测数据N(ωk)最接近。因此关键是构造准则函数，通过优化准则函数，推导出估计系统参数的递推算法或迭代算法。
由文献[7]可知,直接构造实频误差Re[G(jωk)]-M(ωk)平方(和)或虚频误差
Im[G(jωk)]-N(ωk)平方(和)准则函数，表达式太复杂。这里为简化，先采用模型变换，再推导迭代辨识方法。由式(6)可得 Re[G(jω)]+jIm[G(jω)]= (7)
两边同乘以[1+ψT(ω)θ2]+jφT(ω)θ1,比较实部和虚部得到 Re[G(jω)][1+ψT(ω)θ2]-
Im[G(jω)]φT(ω)θ1=b0+ψT(ω)ρ2, (8)
Re[G(jω)]φT(ω)θ1+
Im[G(jω)][1+ψT(ω)θ2]=φT(ω)ρ1。 (9)
实部表达式(8)中包含未知参数b0，θ1,θ2和ρ2，但不包含ρ1,而虚部表达式(9)中包含参数θ1,θ2,ρ1,但不包含b0及ρ2。为此可以定义两个准则函数，一个求解未知参数b0,θ1,θ2和ρ2，另一个求解未知参数ρ1，本节讨论这种方法;也可一个准则函数求解未知参数θ1,θ2,ρ1，另一个求解未知参数b0及ρ2，将在第3节讨论。还可以定义一个联合实部与虚部误差平方和准则函数，将在第4节讨论。当然也可以定义两个准则函数，研究实部虚部耦合辨识方法，将在第5节讨论。 2.1 实部虚部联立梯度迭代估计算法
对系统进行频率特性实验，得到输入正弦信号频率为ωk时对应的频率特性观测数据为H(jωk)=M(ωk)+jN(ωk)。由于存在干扰噪声，它不可能与传递函数模型频率响应特性G(jωk)=Re[G(jωk)]+jIm[G(jωk)]完全相等。为此根据实部表达式(8)，定义包含未知参数向量ϑ的实部方程误差 vr(ωk)=b0+ψT(ωk)ρ2-
M(ωk)[1+ψT(ωk)θ2]+N(ωk)φT(ωk)θ1∈, (10)
根据虚部表达式(9)，定义包含未知参数向量ρ1的虚部方程误差为 vi(ωk)=φT(ωk)ρ1-M(ωk)φT(ωk)θ1- N(ωk)[1+ψT(ωk)θ2]∈。 (11)
定义信息向量
φ(ωk)=
[1,N(ωk)φT(ωk),-M(ωk)ψT(ωk),ψT(ωk)]T。 (12)
考虑从第1个输入频率ω1到第L个输入频率ωL的L组数据，定义堆积实频向量ξr(ωL)，堆积虚频向量ξi(ωL)和堆积信息矩阵Ω(ωL)，Φ(ωL)，Φr(ωL)和Ψi(ωL)如下: ξr(ωL)=
[M(ω1),M(ω2),…,M(ωL)]T∈L, (13) ξi(ωL)=
[N(ω1),N(ω2),…,N(ωL)]T∈L, (14)
Ω(ωL)=[φ(ω1),φ(ω2),…,φ(ωL)], (15)
Φ(ωL)=[φ(ω1),φ(ω2),…,φ(ωL)], (16) Φr(ωL)=
[M(ω1)φ(ω1),M(ω2)φ(ω2),…,M(ωL)φ(ωL)], (17) Ψi(ωL)=
[N(ω1)ψ(ω1),N(ω2)ψ(ω2),…,N(ωL)ψ(ωL)]。 (18)
建立包含L组数据的静态准则函数: J1(ϑ
Vr(ωL)=-ξr(ωL)+ΩT(ωL)ϑ∈L, Vi(ωL)=-ξi(ωL)+ L。
令l=1,2,3,…为迭代变量，和是参数向量ϑ和ρ1在第l次迭代得到的参数估计，μ和μ1为迭代步长(收敛因子)。使用负梯度搜索极小化准则函数J1(ϑ)和J2(ρ1)可以得到迭代关系: (19) (20) (21) (22) (23) (24)
式(20)和(22)分别可以看作状态为和的离散时间系统，为了保证参数向量和收敛，矩阵[I-μΩ(ωL)ΩT(ωL)]和[I-μ1Φ(ωL)ΦT(ωL)]的所有特征值必须在单位圆内，即-I≤I-μΩ(ωL)·ΩT(ωL)≤I和-I≤I-μ1Φ(ωL)ΦT(ωL)≤I，或即0≤μΩ(ωL)ΩT(ωL)≤2I和0≤μ1Φ(ωL)ΦT(ωL)≤2I。因此μ和μ1的一个保守选择是满足
或为计算简单，取为
(25) (26)
联立式(19)，(21)，(23)～(26)，和(12)～(18)，就得到估计系统传递函数(1)参数b0，θ1，θ2，ρ1和ρ2的实部虚部联立梯度迭代算法(joint gradient-based iterative algorithm，J-GI算法): (27)
μ≤2‖Ω(ωL)‖-2, (28) (29) (30)
μ1≤2‖Φ(ωL)‖-2, (31) (32) φ(ωk)=
[1,N(ωk)φT(ωk),-M(ωk)ψT(ωk),ψT(ωk)]T, (33)
ξr(ωL)=[M(ω1),M(ω2),…,M(ωL)]T, (34)
ξi(ωL)=[N(ω1),N(ω2),…,N(ωL)]T, (35)
Ω(ωL)=[φ(ω1),φ(ω2),…,φ(ωL)], (36)
Φ(ωL)=[φ(ω1),φ(ω2),…,φ(ωL)], (37) Φr(ωL)=
[M(ω1)φ(ω1),M(ω2)φ(ω2),…,M(ωL)φ(ωL)], (38) Ψi(ωL)=
[N(ω1)ψ(ω1),N(ω2)ψ(ω2),…,N(ωL)ψ(ωL)], (39) (40)
实部虚部联立梯度迭代算法(27)～(40)计算系统传递函数参数估计和的步骤如下。 1)初始化:令l=1。置初值和为任意实向量，设定数据长度L和参数估计精度ε。 2)施加正弦激励信号，收集实频特性观测数据M(ωk)和虚频特性观测数据N(ωk),k=1,2,…,L。用式(34)～(35)构造堆积实频向量ξr(ωL)和堆积虚频向量ξi(ωL)。
3)视n为奇数或偶数，用式(2)～(3)或式(4)～(5)计算信息向量φ(ωk)和ψ(ωk)，用式(33)构造信息向量φ(ωk)，k=1,2,…,L。
4)用式(36)构造堆积信息矩阵Ω(ωL),用式(37)～(39)构造堆积信息矩阵Φ(ωL)，Φr(ωL)和Ψi(ωL)。
5)根据式(28)和式(31)选择尽量大的迭代步长μ和μ1。
6)用式(29)计算新息用式(32)计算新息
7)用式(27)计算参数估计用式(30)计算参数估计
8)如果就增加1，转至步骤6);否则输出参数估计和结束迭代计算过程。 当然，可以在实部虚部联立梯度迭代算法(27)～(40)中引入加权因子(weighting factor)和遗忘因子(forgetting factor)，得到加权J-GI算法和遗忘因子J-GI算法。 2.2 实部虚部联立多新息梯度迭代算法
设正整数p为新息长度。借助于多新息辨识理论，基于实部虚部联立梯度迭代算法(27)～(40),定义堆积实频向量ξr(p,ωk)，堆积虚频向量ξi(p,ωk)和堆积信息矩阵Ω(p,ωk)，Φ(p,ωk)，Φr(p,ωk)和Ψi(p,ωk)如下: ξr(p,ωk)=
[M(ωk),M(ωk-1),…,M(ωk-p+1)]T∈p, (41) ξi(p,ωk)=
[N(ωk),N(ωk-1),…,N(ωk-p+1)]T∈p, (42) Ω(p,ωk)=
[φ(ωk),φ(ωk-1),…,φ(ωk-p+1)], (43) Φ(p,ωk)=
[φ(ωk),φ(ωk-1),…,φ(ωk-p+1)], (44)
Φr(p,ωk)=[M(ωk)φ(ωk),
M(ωk-1)φ(ωk-1),…,M(ωk-p+1)φ(ωk-p+1)], (45)
Ψi(p,ωk)=[N(ωk)ψ(ωk),
N(ωk-1)ψ(ωk-1),…,N(ωk-p+1)ψ(ωk-p+1)]。 (46)
建立包含最近p组数据的滑动数据窗准则函数: J3(ϑ
Vr(p,ωk)=-ξr(p,ωk)+ΩT(p,ωk)ϑ∈p, Vi(p,ωk)=-ξi(p,ωk)+ΦT(p,ωk)ρ1- p。
令l=1,2,3,…为迭代变量，和是参数向量ϑ和ρ1在角频率ωk时的第l次迭代参数估计。按照2.1节的推导，使用负梯度搜索极小化准则函数J3(ϑ)和J4(ρ1),取收敛因子 (47) (48)
可以得到估计系统传递函数(1)参数b0，θ1，θ2，ρ1和ρ2的实部虚部联立多新息梯度迭代算法(joint multi-innovation gradient-based iterative algorithm，J-MIGI算法): (49) Er,l(p,ωk)=
(50) (51)
Ei,l(p,ωk)=-ξi(p,ωk)+ (52) φ(ωk)=
[1,N(ωk)φT(ωk),-M(ωk)ψT(ωk),ψT(ωk)]T, (53) ξr(p,ωk)=
[M(ωk),M(ωk-1),…,M(ωk-p+1)]T, (54) ξi(p,ωk)=
[N(ωk),N(ωk-1),…,N(ωk-p+1)]T, (55) Ω(p,ωk)=
[φ(ωk),φ(ωk-1),…,φ(ωk-p+1)], (56) Φ(p,ωk)=
[φ(ωk),φ(ωk-1),…,φ(ωk-p+1)], (57)
Φr(p,ωk)=[M(ωk)φ(ωk),
M(ωk-1)φ(ωk-1),…,M(ωk-p+1)φ(ωk-p+1)], (58)
Ψi(p,ωk)=[N(ωk)ψ(ωk),
N(ωk-1)ψ(ωk-1),…,N(ωk-p+1)ψ(ωk-p+1)], (59) (60)
实部虚部联立多新息梯度迭代算法(49)～(60)计算参数估计向量和的步骤如下。 1)初始化:设定新息长度p，参数估计精度ε和最大迭代次数lmax。令k=p，收集实频特性观测数据M(ωj)和虚频特性观测数据N(ωj)，视n为奇数或偶数，用式(2)～(3)或式(4)～(5)计算信息向量φ(ωj)和ψ(ωj)，用式(53)构造信息向量φ(ωj)，j=1,2,…,p-1。置初值和为任意实向量。
2)令l=1。收集实频特性观测数据M(ωk)和虚频特性观测数据N(ωk)，用式(54)～(55)构造堆积实频向量ξr(p,ωk)和堆积虚频向量ξi(p,ωk)。
3)视n为奇数或偶数，用式(2)～(3)或式(4)～(5)计算信息向量φ(ωk)和ψ(ωk)，用式(57)～(59)构造堆积信息矩阵Φ(p,ωk),Φr(p,ωk)和Ψi(p,ωk)。 4)用式(53)构造信息向量φ(ωk)，用式(56)构造堆积信息矩阵Ω(p,ωk)。 5)用式(50)计算新息向量Er,l(p,ωk),用式(52)计算新息向量Ei,l(p,ωk)。 6)用式(49)计算参数估计用式(51)计算参数估计
7)如果l [M(ωk),λM(ωk-1),…,λp-1M(ωk-p+1)]T, (61) ξi(p,ωk)= [N(ωk),λN(ωk-1),…,λp-1N(ωk-p+1)]T, (62) Ω(p,ωk)= [φ(ωk),λφ(ωk-1),…,λp-1φ(ωk-p+1)], (63) Φ(p,ωk)= [φ(ωk),λφ(ωk-1),…,λp-1φ(ωk-p+1)], () Φr(p,ωk)= [M(ωk)φ(ωk),λM(ωk-1)φ(ωk-1),…, λp-1M(ωk-p+1)φ(ωk-p+1)], (65) Ψi(p,ωk)= [N(ωk)ψ(ωk),λN(ωk-1)ψ(ωk-1),…, λp-1N(ωk-p+1)ψ(ωk-p+1)]。 (66) 便得到实部虚部联立遗忘因子多新息梯度迭代算法(joint forgetting factor MIGI algorithm，J-FF-MIGI算法)。 2.3 实部虚部联立最小二乘参数估计算法 令和是参数向量ϑ和ρ1的估计。设ϑ和使准则函数J1(ϑ)和J2(ρ1)达到最小，分别令其对参数向量ϑ和ρ1的偏导数为零，联立式(33)～(40)，就得到估计系统传递函数(1)参数b0，θ1，θ2，ρ1和ρ2的实部虚部联立最小二乘算法(joint least squares algorithm，J-LS算法): (67) (68) φ(ωk)=[1,N(ωk)φT(ωk), -M(ωk)ψT(ωk),ψT(ωk)]T, (69) ξr(ωL)=[M(ω1),M(ω2),…,M(ωL)]T, (70) ξi(ωL)=[N(ω1),N(ω2),…,N(ωL)]T, (71) Ω(ωL)=[φ(ω1),φ(ω2),…,φ(ωL)], (72) Φ(ωL)=[φ(ω1),φ(ω2),…,φ(ωL)], (73) Φr(ωL)=[M(ω1)φ(ω1),M(ω2)φ(ω2),…, M(ωL)φ(ωL)], (74) Ψi(ωL)=[N(ω1)ψ(ω1),N(ω2)ψ(ω2),…, N(ωL)ψ(ωL)], (75) (76) 实部虚部联立最小二乘算法(67)～(76)计算系统传递函数参数估计和的步骤如下。 1)初始化:设定数据长度L。 2)施加角频率为ωk的正弦激励信号，收集实频特性观测数据M(ωk)和虚频特性观测数据N(ωk),k=1,2,…,L。用式(70)～(71)构造堆积实频向量ξr(ωL)和堆积虚频向量ξi(ωL)。 3)视n为奇数或偶数，用式(2)～(3)或式(4)～(5)计算信息向量φ(ωk)和ψ(ωk)，用式(69)构造信息向量φ(ωk)，k=1,2,…,L。 4)用式(72)构造堆积信息矩阵Ω(ωL),用式(73)～(75)构造堆积信息矩阵Φ(ωL)，Φr(ωL)和Ψi(ωL)。 5)用式(67)计算参数估计从式(76)的中读出参数估计和用式(68)计算参数估计 6)输出参数估计和计算结束。 2.4 实部虚部联立多新息最小二乘估计算法 设和是参数向量ϑ和ρ1在角频率ωk时的参数估计。设正整数p为新息长度。借助于多新息辨识理论，基于实部虚部联立最小二乘算法(67)～(76),采用数据长度为p的最新滑动数据窗数据定义的准则函数J3(ϑ)和J4(ρ1),按照2.3节的推导，极小化准则函数J3(ϑ)和J4(ρ1)，分别令其对参数向量ϑ和ρ1的偏导数为零，联立式(53)～(59)，可以得到估计系统传递函数(1)参数b0，θ1，θ2，ρ1和ρ2的实部虚部联立多新息最小二乘算法(joint multi-innovation least squares algorithm，J-MILS算法): [Ω(p,ωk)ΩT(p,ωk)]-1Ω(p,ωk)ξr(p,ωk), (77) [Φ(p,ωk)ΦT(p,ωk)]-1Φ(p,ωk)[ξi(p,ωk)+ (78) φ(ωk)= [1,N(ωk)φT(ωk),-M(ωk)ψT(ωk),ψT(ωk)]T, (79) ξr(p,ωk)= [M(ωk),M(ωk-1),…,M(ωk-p+1)]T, (80) ξi(p,ωk)= [N(ωk),N(ωk-1),…,N(ωk-p+1)]T, (81) Ω(p,ωk)= [φ(ωk),φ(ωk-1),…,φ(ωk-p+1)], (82) Φ(p,ωk)= [φ(ωk),φ(ωk-1),…,φ(ωk-p+1)], (83) Φr(p,ωk)=[M(ωk)φ(ωk), M(ωk-1)φ(ωk-1),…,M(ωk-p+1)φ(ωk-p+1)], (84) Ψi(p,ωk)=[N(ωk)ψ(ωk), N(ωk-1)ψ(ωk-1),…,N(ωk-p+1)ψ(ωk-p+1)], (85) (86) 实部虚部联立多新息最小二乘算法(77)～(86)计算参数估计向量和的步骤如下。 1)初始化:设定新息长度p(p≫n)。令k=p，收集实频特性观测数据M(ωj)和虚频特性观测数据N(ωj)，视n为奇数或偶数，用式(2)～(3)或式(4)～(5)计算信息向量φ(ωj)和ψ(ωj)，用式(79)构造信息向量φ(ωj)，j=1,2,…,p-1。 2)收集实频特性观测数据M(ωk)和虚频特性观测数据N(ωk)，用式(80)～(81)构造堆积实频向量ξr(p,ωk)和堆积虚频向量ξi(p,ωk)。 3)视n为奇数或偶数，用式(2)～(3)或式(4)～(5)计算信息向量φ(ωk)和ψ(ωk)，用式(83)～(85)构造堆积信息矩阵Φ(p,ωk)，Φr(p,ωk)和Ψi(p,ωk)。 4)用式(79)构造信息向量φ(ωk)，用式(82)构造堆积信息矩阵Ω(p,ωk)。 5)用式(77)计算参数估计从式(86)的中读出参数估计和用式(78)计算参数估计 6)输出参数估计增加1，转至步骤2)。 3 联立虚部实部迭代辨识方法 3.1 虚部实部联立梯度迭代估计算法 定义未知参数向量: ϑ (87) 定义信息向量: φ(ωk)=[-M(ωk)φT(ωk), -N(ωk)ψT(ωk),φT(ωk)]T, (88) χ(ωk)=[1,ψT(ωk)]T。 () 那么实部方程误差(10)和虚部方程误差(11)可以表示为 vr(ωk)= b0+ψT(ωk)ρ2-M(ωk)[1+ψT(ωk)θ2]+N(ωk)φT(ωk)θ1=χT(ωk)τ+N(ωk)φT(ωk)θ1-M(ωk)[1+ψT(ωk)θ2]∈, (90) vi(ωk)=φT(ωk)ρ1-M(ωk)φT(ωk)θ1- N(ωk)[1+ψT(ωk)θ2]= -N(ωk)+φT(ωk)ϑ∈。 (91) 考虑从第1个输入频率ω1到第L个输入频率ωL的L组数据，定义堆积实频向量ξr(ωL)，堆积虚频向量ξi(ωL),堆积信息矩阵Ω(ωL)，Φi(ωL)，Ψr(ωL)和Ψ1(ωL)如下: ξr(ωL)= [M(ω1),M(ω2),…,M(ωL)]T∈L, (92) ξi(ωL)= [N(ω1),N(ω2),…,N(ωL)]T∈L, (93) Ω(ωL)=[φ(ω1),φ(ω2),…,φ(ωL)], (94) Φi(ωL)= [N(ω1)φ(ω1),N(ω2)φ(ω2),…,N(ωL)φ(ωL)], (95) Ψr(ωL)= [M(ω1)ψ(ω1),M(ω2)ψ(ω2),…,M(ωL)ψ(ωL)], (96) Ψ1(ωL)=[χ(ω1),χ(ω2),…,χ(ωL)]。 (97) 利用L组数据，分别定义包含未知参数向量和的准则函数: J5(ϑ Vi(ωL)=-ξi(ωL)+ΩT(ωL)ϑ, 令l=1,2,3,…为迭代变量，和是参数向量ϑ和在第l次迭代得到的参数估计，μ和μ1为收敛因子。使用负梯度搜索极小化准则函数J5(ϑ)和J6(τ)，联立式(88)～()和(92)～(97)，就得到估计系统传递函数(1)参数b0，θ1，θ2，ρ1和ρ2的虚部实部联立梯度迭代算法(J-GI算法): (98) (99) μ≤2‖Ω(ωL)‖-2, (100) (101) (102) μ1≤2‖Ψ1(ωL)‖-2， (103) φ(ωk)=[-M(ωk)φT(ωk), -N(ωk)ψT(ωk),φT(ωk)]T, (104) χ(ωk)=[1,ψT(ωk)]T, (105) ξr(ωL)=[M(ω1),M(ω2),…,M(ωL)]T, (106) ξi(ωL)=[N(ω1),N(ω2),…,N(ωL)]T, (107) Ω(ωL)=[φ(ω1),φ(ω2),…,φ(ωL)], (108) Φi(ωL)=[N(ω1)φ(ω1),N(ω2)φ(ω2),…, N(ωL)φ(ωL)], (109) Ψr(ωL)=[M(ω1)ψ(ω1),M(ω2)ψ(ω2),…,M(ωL)ψ(ωL)], (110) Ψ1(ωL)=[χ(ω1),χ(ω2),…,χ(ωL)], (111) (112) (113) 虚部实部联立梯度迭代算法(98)～(113)计算系统传递函数参数估计和的步骤如下。 1)初始化:令l=1。置初值和为任意实向量，设定数据长度L和参数估计精度ε。 2)施加正弦激励信号，收集实频特性观测数据M(ωk)和虚频特性观测数据N(ωk),k=1,2,…,L。用式(106)～(107)构造堆积实频向量ξr(ωL)和堆积虚频向量ξi(ωL)。 3)视n为奇数或偶数，用式(2)～(3)或式(4)～(5)计算信息向量φ(ωk)和ψ(ωk)，用式(104)～(105)构造信息向量φ(ωk)和χ(ωk)，k=1,2,…,L。 4)用式(108)构造堆积信息矩阵Ω(ωL),用式(109)～(111)构造堆积信息矩阵Φi(ωL)，Ψr(ωL)和Ψ1(ωL)。 5)根据式(100)和式(103)选择尽量大的迭代步长μ和μ1。 6)用式(99)计算新息用式(102)计算新息 7)用式(98)计算参数估计从式(112)的中读出参数估计和用式(101)计算参数估计从式(113)的中读出参数估计和 8)如果就增加1，转至步骤6);否则输出参数估计和结束迭代计算过程。 3.2 虚部实部联立多新息梯度迭代算法 设正整数p为新息长度。借助于多新息辨识理论，基于虚部实部联立梯度迭代算法(98)～(113),定义堆积实频向量ξr(p,ωk)，堆积虚频向量ξi(p,ωk),和堆积信息矩阵Ω(p,ωk)，Φi(p,ωk)，Ψr(p,ωk)和Ψ1(p,ωk)如下: ξr(p,ωk)=[M(ωk),M(ωk-1),…, M(ωk-p+1)]T∈p, (114) ξi(p,ωk)=[N(ωk),N(ωk-1),…, N(ωk-p+1)]T∈p, (115) Ω(p,ωk)=[φ(ωk),φ(ωk-1),…, φ(ωk-p+1)], (116) Φi(p,ωk)=[N(ωk)φ(ωk), N(ωk-1)φ(ωk-1),…,N(ωk-p+1)φ(ωk-p+1)], (117) Ψr(p,ωk)=[M(ωk)ψ(ωk), M(ωk-1)ψ(ωk-1),…,M(ωk-p+1)ψ(ωk-p+1)], (118) Ψ1(p,ωk)=[χ(ωk),χ(ωk-1),…,χ(ωk-p+1)]。 (119) 利用最近p组数据，定义包含未知参数向量和的准则函数: J7(ϑ Vi(p,ωk)=-ξi(p,ωk)+ΩT(p,ωk)ϑ, 令l=1,2,3,…为迭代变量，和是参数向量ϑ和在角频率ωk时的第l次迭代参数估计,μ(ωk)和μ1(ωk)为收敛因子。使用负梯度搜索极小化准则函数J7(ϑ)和J8(τ)，联立式(114)～(119),便得到虚部实部联立多新息梯度迭代算法(J-MIGI算法): (120) Ei,l(p,ωk)= (121) μ(ωk)≤2‖Ω(p,ωk)‖-2, (122) (123) Er,l(p,ωk)=-ξr(p,ωk)+ (124) μ1(ωk)≤2‖Ψ1(p,ωk)‖-2, (125) φ(ωk)=[-M(ωk)φT(ωk), -N(ωk)ψT(ωk),φT(ωk)]T, (126) χ(ωk)=[1,ψT(ωk)]T, (127) ξr(p,ωk)=[M(ωk),M(ωk-1),…, M(ωk-p+1)]T, (128) ξi(p,ωk)=[N(ωk),N(ωk-1),…, N(ωk-p+1)]T, (129) Ω(p,ωk)=[φ(ωk),φ(ωk-1),…, φ(ωk-p+1)], (130) Φi(p,ωk)=[N(ωk)φ(ωk), N(ωk-1)φ(ωk-1),…,N(ωk-p+1)φ(ωk-p+1)], (131) Ψr(p,ωk)=[M(ωk)ψ(ωk), M(ωk-1)ψ(ωk-1),…,M(ωk-p+1)ψ(ωk-p+1)], (132) Ψ1(p,ωk)=[χ(ωk),χ(ωk-1),…,χ(ωk-p+1)], (133) (134) (135) 虚部实部联立多新息梯度迭代算法(120)～(135)计算参数估计向量和步骤如下。 1)初始化:设定新息长度p，参数估计精度ε和最大迭代次数lmax。令k=p，收集实频特性观测数据M(ωj)和虚频特性观测数据N(ωj)，视n为奇数或偶数，用式(2)～(3)或式(4)～(5)计算信息向量φ(ωj)和ψ(ωj)，用式(126)～(127)构造信息向量φ(ωj)和χ(ωj)，j=1,2,…,p-1。置初值和为任意实向量。 2)令l=1。收集实频特性观测数据M(ωk)和虚频特性观测数据N(ωk)，用式(128)～(129)构造堆积实频向量ξr(p,ωk)和堆积虚频向量ξi(p,ωk)。 3)视n为奇数或偶数，用式(2)～(3)或式(4)～(5)计算信息向量φ(ωk)和ψ(ωk)，用式(131)～(132)构造堆积信息矩阵Φi(p,ωk)和Ψr(p,ωk)。 4)用式(126)～(127)构造信息向量φ(ωk)和χ(ωk)，用式(130)和(133)构造堆积信息矩阵Ω(p,ωk)和Ψ1(p,ωk)。根据式(122)和(125)确定尽量大的收敛因子μ(ωk)和μ1(ωk)。 5)用式(121)计算新息向量Ei,l(p,ωk),用式(124)计算新息向量Er,l(p,ωk)。 6)用式(120)计算参数估计从式(134)的中读出参数估计和用式(123)计算参数估计从式(135)的中读出参数估计和 7)如果l 令和是参数向量ϑ和的最小二乘估计。极小化准则函数J5(ϑ)和J6(τ)，分别令其对参数向量ϑ和τ的偏导数为零，联立式(104)～(113)，就得到估计系统传递函数(1)参数b0，θ1，θ2，ρ1和ρ2的虚部实部联立最小二乘算法(J-LS算法): (136) (137) φ(ωk)=[-M(ωk)φT(ωk), -N(ωk)ψT(ωk),φT(ωk)]T, (138) χ(ωk)=[1,ψT(ωk)]T, (139) ξr(ωL)=[M(ω1),M(ω2),…,M(ωL)]T, (140) ξi(ωL)=[N(ω1),N(ω2),…,N(ωL)]T, (141) Ω(ωL)=[φ(ω1),φ(ω2),…,φ(ωL)], (142) Φi(ωL)=[N(ω1)φ(ω1),N(ω2)φ(ω2),…, N(ωL)φ(ωL)], (143) Ψr(ωL)=[M(ω1)ψ(ω1),M(ω2)ψ(ω2),…,M(ωL)ψ(ωL)], (144) Ψ1(ωL)=[χ(ω1),χ(ω2),…,χ(ωL)], (145) (146) 虚部实部联立最小二乘算法(136)～(146)计算系统传递函数参数估计和的步骤如下。 1)初始化:设定数据长度L。 2)施加角频率ωk的正弦激励信号，收集实频特性观测数据M(ωk)和虚频特性观测数据N(ωk),k=1,2,…,L。用式(140)～(141)构造堆积实频向量ξr(ωL)和堆积虚频向量ξi(ωL)。 3)视n为奇数或偶数，用式(2)～(3)或式(4)～(5)计算信息向量φ(ωk)和ψ(ωk)，用式(138)～(139)构造信息向量φ(ωk)和χ(ωk)，k=1,2,…,L。 4)用式(142)构造堆积信息矩阵Ω(ωL),用式(143)～(145)构造堆积信息矩阵Φi(ωL)，Ψr(ωL)和Ψ1(ωL)。 5)用式(136)计算参数估计从式(146)的中读出参数估计和用式(137)计算参数估计从式(146)的中读出参数估计和 6)输出参数估计和计算结束。 3.4 虚部实部联立多新息最小二乘估计算法 设正整数p为新息长度。设和是参数向量ϑ和在角频率ωk时的参数估计。借助于多新息辨识理论，基于虚部实部联立最小二乘算法(136)～(146),采用数据长度为p的最新滑动数据窗数据定义的准则函数J7(ϑ)和J8(τ),按照3.3节的推导，极小化准则函数J7(ϑ)和J8(τ)，分别令其对参数向量ϑ和τ的偏导数为零，联立式(126)～(133)，便得到实部虚部联立多新息最小二乘算法(J-MILS算法): [Ω(p,ωk)ΩT(p,ωk)]-1Ω(p,ωk)ξi(p,ωk), (147) (148) φ(ωk)=[-M(ωk)φT(ωk), -N(ωk)ψT(ωk),φT(ωk)]T, (149) χ(ωk)=[1,ψT(ωk)]T, (150) ξr(p,ωk)=[M(ωk),M(ωk-1),…, M(ωk-p+1)]T, (151) ξi(p,ωk)=[N(ωk),N(ωk-1),…, N(ωk-p+1)]T, (152) Ω(p,ωk)=[φ(ωk),φ(ωk-1),…, φ(ωk-p+1)], (153) Φi(p,ωk)= [N(ωk)φ(ωk),N(ωk-1)φ(ωk-1),…, N(ωk-p+1)φ(ωk-p+1)], (154) Ψr(p,ωk)=[M(ωk)ψ(ωk), M(ωk-1)ψ(ωk-1),…,M(ωk-p+1)ψ(ωk-p+1)], (155) Ψ1(p,ωk)=[χ(ωk),χ(ωk-1),…, χ(ωk-p+1)], (156) (157) (158) 虚部实部联立多新息最小二乘算法(147)～(158)计算参数估计向量和的步骤如下。 1)初始化:设定新息长度p(p≫n)。令k=p，收集实频特性观测数据M(ωj)和虚频特性观测数据N(ωj)，视n为奇数或偶数，用式(2)～(3)或式(4)～(5)计算信息向量φ(ωj)和ψ(ωj)，用式(149)～(150)构造信息向量φ(ωj)和χ(ωj)，j=1,2,…,p-1。 2)收集实频特性观测数据M(ωk)和虚频特性观测数据N(ωk)，用式(151)～(152)构造堆积实频向量ξr(p,ωk)和堆积虚频向量ξi(p,ωk)。 3)视n为奇数或偶数，用式(2)～(3)或式(4)～(5)计算信息向量φ(ωk)和ψ(ωk)，用式(149)～(150)构造信息向量φ(ωk)和χ(ωk)。 4)用式(153)构造堆积信息矩阵Ω(p,ωk)，用式(154)～(156)构造堆积信息矩阵Φi(p,ωk)，Ψr(p,ωk)和Ψ1(p,ωk)。 5)用式(147)计算参数估计从式(157)的中读出参数估计和用式(148)计算参数估计从式(158)的中读出参数估计和 6)输出参数估计和增加1，转至步骤2)。 4 实部虚部联合迭代辨识方法 实部方程误差(10)中包含未知参数b0，θ1,θ2和ρ2，但不包含ρ1,虚部方程误差 式(11)中包含参数θ1,θ2，ρ1，但不包含b0及ρ2。因此无法仅凭借其中的一个关系式同时获得系统全部参数的估计。为此需要定义一个联合实部方程误差与虚部方程误差平方和准则函数。本节研究实部虚部联合梯度迭代算法、实部虚部联合多新息梯度迭代算法、实部虚部联合最小二乘算法、实部虚部联合多新息最小二乘算法。 4.1 联合梯度迭代估计算法 为了估计出系统的参数b0,θ1,θ2,ρ1和ρ2，定义参数向量ϑ根据式(10)～(11),构造包含误差vr(ωk)和vi(ωk)的平方和准则函数 J9(ϑ)=J9(b0,θ1,θ2,ρ1,ρ2)= 其中L为数据长度。定义实频虚频观测数据向量ξ(ωk)和信息矩阵γ(ωk)如下: ξ(ωk)=[-M(ωk),-N(ωk)]T∈2, (159) γ(ωk)= (160) 考虑从第1个输入频率ω1到第L个输入频率ωL的L组数据，定义联合堆积向量Ξ(ωL)和堆积信息矩阵Γ(ωL)如下: 2L, (161) Γ(ωL)=[γT(ω1),γT(ω2),…,γT(ωL)]。 (162) 因此有 J9(ϑϑ‖2。 令l=1,2,3,…为迭代变量，表示系统参数ϑ的第l次迭代估计，μl为迭代步长(收敛因子)。使用负梯度搜索，极小化准则函数J9(ϑ)得到迭代关系: (163) 其中 (1) 下面求最佳步长。将ϑ代入准则函数J9(ϑ)得到 极小化g(μl),令得到 -2‖Γ(ωL)Εl(ωL)‖2+ 2μl‖ΓT(ωL)Γ(ωL)Εl(ωL)‖2=0。 由此，给出最佳收敛因子 (165) 进一步简化，收敛因子取为 (166) 联立式(159)～(166),可以得到估计系统传递函数(1)参数b0,θ1,θ2,ρ1和ρ2的实部虚部联合梯度迭代算法(J-GI算法): (167) (168) (169) 或 μl=‖Γ(ωL)‖-2, (170) ξ(ωk)=[-M(ωk),-N(ωk)]T, (171) Υ(ωk)= (172) Ξ(ωL)=[ξT(ω1),ξT(ω2),…,ξT(ωL)]T, (173) Γ(ωL)=[γT(ω1),γT(ω2),…,γT(ωL)], (174) (175) 实部虚部联合梯度迭代算法(167)～(175)计算系统传递函数参数估计的步骤如下。 1)初始化:令l=1。置初值为任意实向量，设定数据长度L(L≫n)和参数估计精度ε。 2)施加正弦激励信号，收集实频特性观测数据M(ωk)和虚频特性观测数据N(ωk)，用式(171)构成联合向量ξ(ωk),k=1,2,…,L。用式(173)构造联合堆积实频虚频向量Ξ(ωL)。 3)视n为奇数或偶数，用式(2)～(3)或式(4)～(5)计算信息向量φ(ωk)和ψ(ωk)，k=1,2,…,L。用式(172)构造信息矩阵γ(ωk)。用式(174)构造堆积信息矩阵Γ(ωL)。 4)用式(168)计算新息向量Εl(ωL)。用式(169)或(170)确定收敛因子μl。 5)用式(167)计算参数估计从式(175)的中读出参数估计和 6)如果就增加1，转至步骤4);否则输出参数估计结束迭代计算过程。 4.2 联合多新息梯度迭代算法 设正整数p为新息长度。定义移动数据窗准则函数 J10(ϑ)=J10(b0,θ1,θ2,ρ1,ρ2)= 定义堆积实频虚频观测数据向量Ξ(p,ωk)和堆积信息矩阵Γ(p,ωk)如下: 2p, Γ(p,ωk)=[γT(ωk),γT(ωk-1),…, γT(ωk-p+1)]。 令l=1,2,3,…为迭代变量，表示系统参数ϑ在角频率ωk时的第l次迭代估计，μl(ωk)为在角频率ωk时的收敛因子。借助于多新息辨识理论，参照实部虚部联合梯度迭代算法(167)～(175)的推导,使用负梯度搜索，极小化准则函数J10(ϑ)，可以得到估计系统传递函数(1)参数b0,θ1,θ2,ρ1和ρ2的实部虚部联合多新息梯度迭代算法(J-MIGI算法): (176) (177) (178) 或 μl(ωk)=‖Γ(p,ωk)‖-2, (179) ξ(ωk)=[-M(ωk),-N(ωk)]T, (180) γ(ωk)= (181) Ξ(p,ωk)=[ξT(ωk),ξT(ωk-1),…, ξT(ωk-p+1)]T, (182) Γ(p,ωk)=[γT(ωk),γT(ωk-1),…, γT(ωk-p+1)], (183) (184) 实部虚部联合多新息梯度迭代算法(176)～(184)计算参数估计向量的步骤如下。 1)初始化:设定新息长度p(p≫n)，参数估计精度ε和最大迭代次数lmax。令k=p，置参数估计初值收集实频特性观测数据M(ωj)和虚频特性观测数据N(ωj)，视n为奇数或偶数，用式(2)～(3)或式(4)～(5)计算信息向量φ(ωj)和ψ(ωj)，用式(180)～(181)构造联合实频虚频信息向量ξ(ωj)和信息矩阵γ(ωj)，j=1,2,…,p-1。 2)令l=1,收集实频特性观测数据M(ωk)和虚频特性观测数据N(ωk)。视n为奇数或偶数，用式(2)～(3)或式(4)～(5)计算信息向量φ(ωk)和ψ(ωk)，用式(180)～(181)构造联合实频虚频信息向量ξ(ωk)和信息矩阵γ(ωk)。用式(182)～(183)构 造堆积联合实频虚频向量Ξ(p,ωk)和堆积信息矩阵Γ(p,ωk)。 3)用式(177)计算新息向量Εl(p,ωk)。用式(178)或(179)确定收敛因子μl(ωk)。 4)用式(176)计算参数估计从式(184)的中读出参数估计和 5)如果l 表示统参数ϑ的最小二乘参数估计。极小化准则函数J9(ϑ)，令其对参数向量ϑ的偏导数为零，联立式(171)～(175)，就得到估计系统传递函数(1)参数b0,θ1,θ2,ρ1和ρ2的联合最小二乘算法(J-LS算法): (185) ξ(ωk)=[-M(ωk),-N(ωk)]T, (186) γ(ωk)= (187) Ξ(ωL)=[ξT(ω1),ξT(ω2),…,ξT(ωL)]T, (188) Γ(ωL)=[γT(ω1),γT(ω2),…,γT(ωL)], (1) (190) 实部虚部联合最小二乘算法(185)～(190)计算系统传递函数参数估计的步骤如下。 1)初始化:设定数据长度L(L≫n)。 2)施加正弦激励信号，收集实频特性观测数据M(ωk)和虚频特性观测数据N(ωk)，用式(186)构成联合向量ξ(ωk),k=1,2,…,L。用式(188)构造联合堆积实频虚频向量Ξ(ωL)。 3)视n为奇数或偶数，用式(2)～(3)或式(4)～(5)计算信息向量φ(ωk)和 ψ(ωk),k=1,2,…,L。用式(187)构造信息矩阵γ(ωk)。用式(1)构造堆积信息矩阵Γ(ωL)。 4)用式(185)计算参数估计从式(190)的中读出参数估计和结束计算过程。 4.4 联合多新息最小二乘估计算法 令表示系统参数ϑ在角频率ωk时的估计。设正整数p为新息长度。借助于多新息辨识理论，基于虚部实部联合最小二乘算法(185)～(190),采用数据长度为p的最新滑动数据窗数据定义的准则函数J10(ϑ),按照4.3节的推导，极小化准则函数J10(ϑ)，令其对参数向量ϑ的偏导数为零,联立式(180)～(183)，就得到估计系统传递函数(1)参数b0,θ1,θ2,ρ1和ρ2的实部虚部联合多新息最小二乘算法(J-MILS算法): [Γ(p,ωk)ΓT(p,ωk)]-1Γ(p,ωk)Ξ(p,ωk), (191) ξ(ωk)=[-M(ωk),-N(ωk)]T, (192) γ(ωk)= (193) Ξ(p,ωk)=[ξT(ωk),ξT(ωk-1),…, ξT(ωk-p+1)]T, (194) Γ(p,ωk)=[γT(ωk),γT(ωk-1),…, γT(ωk-p+1)], (195) (196) 实部虚部联合多新息最小二乘算法(191)～(196)计算参数估计向量的步骤如下。 1)初始化:设定新息长度p(p≫n)，设定参数估计精度ε。令k=p，收集实频特性观测数据M(ωj)和虚频特性观测数据N(ωj)，视n为奇数或偶数，用式(2)～(3)或式(4)～(5)计算信息向量φ(ωj)和ψ(ωj)，用式(192)～(193)构造联合实频虚频信息向量ξ(ωj)和信息矩阵Υ(ωj),j=1,2,…,p-1。 2)收集实频特性观测数据M(ωk)和虚频特性观测数据N(ωk)。视n为奇数或偶数，用式(2)～(3)或式(4)～(5)计算信息向量φ(ωk)和ψ(ωk)，用式(192)～(193)构造联合实频虚频信息向量ξ(ωk)和信息矩阵γ(ωk)。用式(194)～(195)构造堆积联合实频虚频向量Ξ(p,ωk)和堆积信息矩阵Γ(p,ωk)。 3)用式(191)计算参数估计从式(196)的中读出参数估计和 4)如果就增加1，转至步骤2);否则输出参数估计结束计算过程。 5 实部虚部耦合迭代辨识方法 定义有关参数向量和信息向量: ϑ (197) φ1(ωk)= [N(ωk)φT(ωk),-M(ωk)ψT(ωk)]T, (198) φ2(ωk)= [-M(ωk)φT(ωk),-N(ωk)ψT(ωk)]T, (199) χ(ωk)=[1,ψT(ωk)]T。 (200) 根据实部误差和虚部误差定义式(10)～(11),我们可以发现:实部误差vr(ωk)和虚部误差vi(ωk)包含了不同的参数向量τ和ρ1，还包含了相同的参数向量ϑ，因此可以用耦合辨识概念来研究基于频率响应的耦合辨识方法。 5.1 耦合梯度迭代估计算法 考虑从第1个输入频率ω1到第L个输入频率ωL的L组数据，定义堆积实频向量ξr(ωL)，堆积虚频向量ξi(ωL)，堆积信息矩阵Ω1(ωL),Ω2(ωL),Φ(ωL)和Ψ1(ωL)如下: ξr(ωL)=[M(ω1),M(ω2),…,M(ωL)]T∈L, (201) ξi(ωL)=[N(ω1),N(ω2),…,N(ωL)]T∈L, (202) Ω1(ωL)=[φ1(ω1),φ1(ω2),…,φ1(ωL)], (203) Ω2(ωL)=[φ2(ω1),φ2(ω2),…,φ2(ωL)], (204) Φ(ωL)=[φ(ω1),φ(ω2),…,φ(ωL)], (205) Ψ1(ωL)=[χ(ω1),χ(ω2),…,χ(ωL)]。 (206) 根据式(10)，定义包含未知参数向量ϑ和τ的准则函数 J11(ϑϑ+ 令l=1,2,3,…为迭代变量,和分别是参数向量ϑ和的第l次迭代估计。极小化准则函数J11(ϑ,τ)，取收敛因子为 (207) (208) 联立式(207)～(208),(198),(200),(201),(203)和(206),就得到估计参数向量ϑ和τ的梯度迭代算法: (209) (210) φ1(ωk)=[N(ωk)φT(ωk),-M(ωk)ψT(ωk)]T, (211) χ(ωk)=[1,ψT(ωk)]T, (212) ξr(ωL)=[M(ω1),M(ω2),…,M(ωL)]T, (213) Ω1(ωL)=[φ1(ω1),φ1(ω2),…,φ1(ωL)], (214) Ψ1(ωL)=[χ(ω1),χ(ω2),…,χ(ωL)], (215) (216) 根据式(11)，定义包含未知参数向量ϑ和ρ1的准则函数 J12(ϑϑ+ ΦT(ωL)ρ1‖2。 令和分别是参数向量ϑ和ρ1的第l次迭代估计。极小化准则函数J12(ϑ,ρ1)，可以得到估计参数向量ϑ和ρ1的梯度迭代算法: (217) (218) φ2(ωk)=[-M(ωk)φT(ωk), -N(ωk)ψT(ωk)]T, (219) ξi(ωL)=[N(ω1),N(ω2),…,N(ωL)]T, (220) Ω2(ωL)=[φ2(ω1),φ2(ω2),…,φ2(ωL)], (221) Φ(ωL)=[φ(ω1),φ(ω2),…,φ(ωL)], (222) (223) 算法(209)～(216)可以估计出ϑ和τ,算法(217)～(223)可以估计出ϑ和ρ1,所以ϑ 被估计了2次,得到了2个估计和因此可以用它们的平均值作为ϑ的估计 (224) 于是联立式(209)～(216)和式(217)～(224),便得到估计参数向量ϑ,τ和ρ1的基于平均的耦合梯度迭代算法: (225) (226) (227) (228) φ1(ωk)=[N(ωk)φT(ωk), -M(ωk)ψT(ωk)]T, (229) φ2(ωk)=[-M(ωk)φT(ωk), -N(ωk)ψT(ωk)]T, (230) χ(ωk)=[1,ψT(ωk)]T, (231) ξr(ωL)=[M(ω1),M(ω2),…,M(ωL)]T, (232) ξi(ωL)=[N(ω1),N(ω2),…,N(ωL)]T, (233) Ω1(ωL)=[φ1(ω1),φ1(ω2),…,φ1(ωL)], (234) Ω2(ωL)=[φ2(ω1),φ2(ω2),…,φ2(ωL)], (235) Φ(ωL)=[φ(ω1),φ(ω2),…,φ(ωL)], (236) Ψ1(ωL)=[χ(ω1),χ(ω2),…,χ(ωL)], (237) (238) (239) 由于两个准则函数J11(ϑ,τ)和J12(ϑ,ρ1)包含了相同的参数向量ϑ，所以利用耦合辨识概念,将基于平均的耦合梯度迭代算法(225)～(239)的式(225)～(228)和(238)～(239)修改为 (240) (241) (242) (243) (244) 就得到估计参数向量ϑ,τ和ρ1的耦合梯度迭代算法(coupled gradient-based iterative algorithm，C-GI算法)。 耦合梯度迭代算法(240)～(244)是先计算和后计算和反过来，先计算和后计算和则只需将耦合梯度迭代算法式(240)～(244)修改为 (245) (246) (247) (248) (249) 5.2 耦合多新息梯度迭代算法 设正整数p为新息长度。借助于多新息辨识理论,定义堆积实频向量ξr(p,ωk)和堆积虚频向量ξi(p,ωk)如下: ξr(p,ωk)=[M(ωk),M(ωk-1),…, M(ωk-p+1)]T∈p, (250) ξi(p,ωk)=[N(ωk),N(ωk-1),…, N(ωk-p+1)]T∈p。 (251) 将信息向量φ1(ωk),φ2(ωk),φ(ωk)和χ(ωk)扩展为堆积信息矩阵Ω1(p,ωk),Ω2(p,ωk),Φ(p,ωk)和Ψ1(p,ωk)如下: Ω1(p,ωk)=[φ1(ωk),φ1(ωk-1),…, φ1(ωk-p+1)], (252) Ω2(p,ωk)=[φ2(ωk),φ2(ωk-1),…, φ2(ωk-p+1)], (253) Φ(p,ωk)=[φ(ωk),φ(ωk-1),…, φ(ωk-p+1)], (254) Ψ1(p,ωk)=[χ(ωk),χ(ωk-1),…, χ(ωk-p+1)]。 (255) 根据式(10)，定义包含未知参数向量ϑ和τ的准则函数 J13(ϑ,τ)= ϑ 令l=1,2,3,…为迭代变量。设和分别是参数向量ϑ和在角频率ωk时的第l次迭代估计,也就是用梯度搜索极小化准则函数J13(ϑ,τ)得到的估计;根据式(11)，定义包含未知参数向量ϑ和ρ1的准则函数 J14(ϑϑ+ΦT(p,ωk)ρ1‖2。 令和分别是参数向量ϑ和ρ1在角频率ωk时的第l次迭代估计,也就是用梯度搜索极小化准则函数J14(ϑ,ρ1)得到的估计;按照5.1节的推导方法，利用耦合辨识概念， 可以得到估计参数向量ϑ,τ和ρ1的耦合多新息梯度迭代算法(coupled multi-innovation gradient-based iterative algorithm，C-MIGI算法): (256) (257) (258) (259) φ1(ωk)= [N(ωk)φT(ωk),-M(ωk)ψT(ωk)]T, (260) φ2(ωk)= [-M(ωk)φT(ωk),-N(ωk)ψT(ωk)]T, (261) χ(ωk)=[1,ψT(ωk)]T, (262) ξr(p,ωk)=[M(ωk),M(ωk-1),…, M(ωk-p+1)]T, (263) ξi(p,ωk)=[N(ωk),N(ωk-1),…, N(ωk-p+1)]T, (2) Ω1(p,ωk)=[φ1(ωk),φ1(ωk-1),…, φ1(ωk-p+1)], (265) Ω2(p,ωk)=[φ2(ωk),φ2(ωk-1),…, φ2(ωk-p+1)], (266) Ω(p,ωk)=[φ(ωk),φ(ωk-1),…,φ(ωk-p+1)], (267) Ψ1(p,ωk)=[χ(ωk),χ(ωk-1),…, χ(ωk-p+1)], (268) (269) (270) 实部虚部耦合多新息梯度迭代算法(256)～(270)计算参数估计向量和的步骤如下。 1)初始化:设定新息长度p(p≫n)，参数估计精度ε和最大迭代次数lmax。令k=p，置参数估计初值和为任意实向量。收集实频特性观测数据M(ωj)和虚频特性观测数据N(ωj)，视n为奇数或偶数，用式(2)～(3)或式(4)～(5)计算信息向量φ(ωj)和ψ(ωj)，用式(260)～(262)构造信息向量φ1(ωj),φ2(ωj)和χ(ωj),j=1,2,…,p-1。 2)令l=1。收集实频特性观测数据M(ωk)和虚频特性观测数据N(ωk)，视n为奇数或偶数，用式(2)～(3)或式(4)～(5)计算信息向量φ(ωk)和ψ(ωk)，用式(260)～(262)构造信息向量φ1(ωk)，φ2(ωk)和χ(ωk)。 3)用式(263)～(2)构造堆积实频向量ξr(p,ωk)和堆积虚频向量ξi(p,ωk),用式 (265)～(268)构造堆积信息矩阵Ω1(p,ωk),Ω2(p,ωk),Φ(p,ωk)和Ψ1(p,ωk)。 4)用式(256)～(257)计算参数估计和用式(258)和(259)计算参数估计和 5)根据式(269)置并从中读取参数估计和从式(270)的中读取参数估计和 6)如果l 令l=1,2,3,…为迭代变量，和分别是参数向量ϑ和的第l次迭代估计。极小化准则函数J11(ϑ,τ)，分别令其对参数向量ϑ和τ的偏导数为零,联立式(211)～(216)，就得到估计参数向量ϑ和τ的最小二乘迭代估计算法: (271) (272) φ1(ωk)= [N(ωk)φT(ωk),-M(ωk)ψT(ωk)]T, (273) χ(ωk)=[1,ψT(ωk)]T, (274) ξr(ωL)=[M(ω1),M(ω2),…,M(ωL)]T, (275) Ω1(ωL)=[φ1(ω1),φ1(ω2),…,φ1(ωL)], (276) Ψ1(ωL)=[χ(ω1),χ(ω2),…,χ(ωL)], (277) (278) 令和分别是参数向量ϑ和ρ1的第l次迭代估计。极小化准则函数J12(ϑ,ρ1)，可以得到估计参数向量ϑ和ρ1的最小二乘迭代估计算法: (279) (280) φ2(ωk)=[-M(ωk)φT(ωk), -N(ωk)ψT(ωk)]T, (281) ξi(ωL)=[N(ω1),N(ω2),…,N(ωL)]T, (282) Ω2(ωL)=[φ2(ω1),φ2(ω2),…,φ2(ωL)], (283) Φ(ωL)=[φ(ω1),φ(ω2),…,φ(ωL)], (284) (285) 算法(271)～(278)可以估计出ϑ和τ,算法(279)～(285)可以估计出ϑ和ρ1,所以ϑ被估计了两次,得到了两个估计和因此可以用它们的平均值作为ϑ的估计 (286) 于是式(271)～(286)构成了估计参数向量ϑ,τ和ρ1的基于平均的耦合最小二乘迭代算法: (287) (288) (2) (290) φ1(ωk)=[N(ωk)φT(ωk), -M(ωk)ψT(ωk)]T, (291) φ2(ωk)=[-M(ωk)φT(ωk), -N(ωk)ψT(ωk)]T, (292) χ(ωk)=[1,ψT(ωk)]T, (293) ξr(ωL)=[M(ω1),M(ω2),…,M(ωL)]T, (294) ξi(ωL)=[N(ω1),N(ω2),…,N(ωL)]T, (295) Ω1(ωL)=[φ1(ω1),φ1(ω2),…,φ1(ωL)], (296) Ω2(ωL)=[φ2(ω1),φ2(ω2),…,φ2(ωL)], (297) Φ(ωL)=[φ(ω1),φ(ω2),…,φ(ωL)], (298) Ψ1(ωL)=[χ(ω1),χ(ω2),…,χ(ωL)], (299) (300) (301) 由于两个准则函数J11(ϑ,τ)和J12(ϑ,ρ1)包含了相同的参数向量ϑ，所以利用耦合辨识概念,将基于平均的耦合梯度迭代算法(287)～(301)的式(287)～(290)和(300)～(301)修改为 (302) (303) (304) (305) (306) 就得到估计参数向量ϑ,τ和ρ1的耦合最小二乘迭代算法(coupled least squares based iterative algorithm，C-LSI算法)。 5.4 耦合多新息最小二乘迭代算法 设正整数p为新息长度。采用数据长度为p的最新滑动数据窗数据定义的准则函数J13(ϑ,τ)和J14(ϑ,ρ1),令l=1,2,3,…为迭代变量。设和分别是参数向量ϑ和在角频率ωk时的第l次迭代估计,也就是极小化准则函数J13(ϑ,τ)，令其对参数向量ϑ和τ的梯度为零得到的最小二乘估计;令和分别是参数向量ϑ和ρ1在角频率ωk时的第l次迭代估计,也就是极小化准则函数J14(ϑ,ρ1)，令其对参数向量ϑ和ρ1的梯度为零,得到的最小二乘估计。 借助于多新息辨识理论和耦合辨识概念，按照5.3节的推导方法，极小化准则函数J13(ϑ,τ)和J14(ϑ,ρ1),联立式(260)～(270)，可以得到估计参数向量ϑ,τ和ρ1的实部虚部耦合多新息最小二乘迭代算法(coupled multi-innovation least squares based iterative algorithm，C-MILSI算法): (307) (308) (309) (310) φ1(ωk)=[N(ωk)φT(ωk), -M(ωk)ψT(ωk)]T, (311) φ2(ωk)=[-M(ωk)φT(ωk), -N(ωk)ψT(ωk)]T, (312) χ(ωk)=[1,ψT(ωk)]T, (313) ξr(p,ωk)=[M(ωk),M(ωk-1),…, M(ωk-p+1)]T, (314) ξi(p,ωk)=[N(ωk),N(ωk-1),…, N(ωk-p+1)]T, (315) Ω1(p,ωk)=[φ1(ωk),φ1(ωk-1),…, φ1(ωk-p+1)], (316) Ω2(p,ωk)=[φ2(ωk),φ2(ωk-1),…, φ2(ωk-p+1)], (317) Φ(p,ωk)=[φ(ωk),φ(ωk-1),…, φ(ωk-p+1)], (318) Ψ1(p,ωk)=[χ(ωk),χ(ωk-1),…, χ(ωk-p+1)], (319) (320) (321) 实部虚部耦合多新息最小二乘迭代算法(307)～(321)计算参数估计向量和的步骤如下。 1)初始化:设定新息长度p(p≫n)，参数估计精度ε和最大迭代次数lmax。令k=p，置参数估计初值和为任意实向量。收集实频特性观测数据M(ωj)和虚频特性观测数据N(ωj)，视n为奇数或偶数，用式(2)～(3)或式(4)～(5)计算信息向量φ(ωj)和ψ(ωj)，用式(311)～(313)构造信息向量φ1(ωj),φ2(ωj)和χ(ωj),j=1,2,…,p-1。 2)令l=1。收集实频特性观测数据M(ωk)和虚频特性观测数据N(ωk)，视n为奇数或偶数，用式(2)～(3)或式(4)～(5)计算信息向量φ(ωk)和ψ(ωk)，用式(311)～(313)构造信息向量φ1(ωk),φ2(ωk)和χ(ωk)。 3)用式(314)～(315)构造堆积实频向量ξr(p,ωk)和堆积虚频向量ξi(p,ωk),用式(316)～(319)构造堆积信息矩阵Ω1(p,ωk),Ω2(p,ωk),Φ(p,ωk)和Ψ1(p,ωk)。 4)用式(307)～(308)计算参数估计和用式(309)～(310)计算参数估计和 5)根据式(320)置并从中读取参数估计和从式(320)的中读取参数估计和 6)如果l (324) (325) (326) (327) (328) 得到基于平均的耦合多新息最小二乘迭代算法。 6 结 语 针对阶次已知的一般线性时不变系统传递函数，利用系统的频率特性观测数据,即实频特性数据和虚频特性数据，根据系统频率特性参数化辨识模型,为避免复杂的非线性优化问题，通过模型转化,定义实部误差和虚部误差，把非线性问题转化为二次优化问题,利用梯度搜索、最小二乘搜索、多新息辨识理论和耦合辨识概念,提出了辨识传递函数参数的梯度迭代算法、多新息梯度迭代算法、最小二乘(迭代)算法、多新息最小二乘(迭代)算法、联合迭代辨识算法、耦合迭代辨识算法等。值得指出的是，该方法对一般线性时不变系统参数辨识都有效，对系统是否含有共轭复数极点、重极点等没有。 参 考 文 献 【相关文献】 [1] 丁锋. 现代控制理论[M]. 北京: 清华大学出版社, 2018. DING Feng. Modern Control Theory [M]. Beijing: Tsinghua University Press, 2018. [2] 丁锋. 系统辨识新论[M]. 北京: 科学出版社, 2013. DING Feng. System Identification—New Theory and Methods [M]. Beijing: Science Press, 2013. [3] 丁锋. 系统辨识—迭代搜索原理与辨识方法[M]. 北京: 科学出版社, 2018. DING Feng. System Identification—Iterative Search Principle and Identification Methods [M]. Beijing: Science Press, 2018. [4] 丁锋. 系统辨识—多新息辨识理论与方法[M]. 北京: 科学出版社, 2016. DING Feng. System Identification—Multi-Innovation Identification Theory and Methods [M]. Beijing: Science Press, 2016. [5] 丁锋, 徐玲, 刘喜梅. 传递函数辨识(3): 正弦响应两点法和多点法[J]. 青岛科技大学学报(自然科学版), 2018, 39(3): 1-15. DING Feng, XU Ling, LIU Ximei. Transfer function identification. Part C: Two-point and multi-point methods based on the sine responses [J].Journal of Qingdao University of Science and Technology (Natural Science Edition), 2018, 39(3): 1-15. [6] 丁锋, 徐玲, 刘喜梅. 传递函数辨识(8): 基于正弦响应的参数估计方法[J]. 青岛科技大学学报(自然科学版), 2018, 40(2): 1-19. DING Feng, XU Ling, LIU Ximei. Transfer function identification. Part H: Parameter estmation methods based on the sine responses [J]. Journal of Qingdao University of Science and Technology (Natural Science Edition), 2018, 40(2): 1-19. [7] 丁锋, 徐玲, 刘喜梅. 传递函数辨识(9): 基于频率响应的递推参数估计方法[J]. 青岛科技大学学报(自然科学版), 2019, 40(4): 1-24. DING Feng, XU Ling, LIU Ximei. Transfer function identification. Part I: Recursive parameter estimation methods based on the frequency responses [J]. Journal of Qingdao University of Science and Technology (Natural Science Edition), 2019, 40(4): 1-24. [8] 陈垒. 连续系统的传递函数模型辨识[D]. 无锡：江南大学, 2012. CHEN Lei. Identification for Transfer Function Models of Continuous-Time Systems [D]. Wuxi: Jiangnan University, 2012. [9] CHEN L, LI J H, DING R. Identification of the second-order systems based on the step response [J]. Mathematical and Computer Modelling, 2011, 53(5/6): 1074-1083. [10] XU L, CHEN L, XIONG W L. Parameter estimation and controller design for dynamic systems from the step responses based on the Newton iteration [J]. Nonlinear Dynamics, 2015, 79(3): 2155-2163. [11] XU L. Application of the Newton iteration algorithm to the parameter estimation for dynamical systems [J]. Journal of Computational and Applied Mathematics, 2015, 288: 33-43. [12] XU L. The damping iterative parameter identification method for dynamical systems based on the sine signal measurement [J]. Signal Processing, 2016, 120: 660-667. [13] XU L, DING F. The parameter estimation algorithms for dynamical response signals based on the multi-innovation theory and the hierarchical principle [J]. IET Signal Processing, 2017, 11(2): 228-237. [14] XU L, DING F. Parameter estimation for control systems based on impulse responses [J]. International Journal of Control Automation and Systems, 2017, 15(6): 2471-2479. [15] 徐玲. 基于移动数据窗的传递函数多新息随机梯度辨识方法[J]. 控制与决策, 2017, 32(6): 1091-1096. XU Ling. Moving data window based multi-innovation stochastic gradient identification method for transfer functions [J]. Control and Decision, 2017, 32(6): 1091-1096. [16] XU L, DING F, ZHU Q M. Hierarchical Newton and least squares iterative estimation algorithm for dynamic systems by transfer functions based on the impulse responses [J]. International Journal of Systems Science, 2019, 50(1): 141-151. [17] XU L. The parameter estimation algorithms based on the dynamical response measurement data [J]. Advances in Mechanical Engineering, 2017, 9(11): 1-12. [18] XU L, DING F. Recursive least squares and multi-innovation stochastic gradient parameter estimation methods for signal modeling [J]. Circuits Systems and Signal Processing, 2017, 36(4): 1735-1753. [19] ZHOU L C, LI X L, XU H G, et al. Multi-innovation stochastic gradient method for harmonic modelling of power signals [J]. IET Signal Processing, 2016, 10(7): 737-742. [20] CAO Y N, LIU Z Q. Signal frequency and parameter estimation for power systems using the hierarchical identification principle [J]. Mathematical and Computer Modelling, 2010, 52(5/6): 854-861. [21] XU L, XIONG W L, ALSAEDI A, et al. Hierarchical parameter estimation for the frequency response based on the dynamical window data [J]. International Journal of Control Automation and Systems, 2018, 16(4): 1756-17. [22] PAN J, JIANG X, WAN X K, et al. A filtering based multi-innovation extended stochastic gradient algorithm for multivariable control systems [J]. International Journal of Control Automation and Systems, 2017, 15(3): 11-1197.
因篇幅问题不能全部显示，请点此查看更多更全内容
Copyright © 2019- 517ttc.cn 版权所有 赣ICP备2024042791号-8
违法及侵权请联系：TEL:199 18 7713 E-MAIL:2724546146@qq.com
本站由北京市万商天勤律师事务所王兴未律师提供法律服务