CN103929151A - 自适应最优相角陷波滤波器设计方法 - Google Patents

Abstract

本发明具体公开了一种自适应最优相角陷波滤波器的设计方法，用于抑制未知频段的窄频带信号干扰，使系统具有抵消低频段的宽频带随机干扰信号特性。其系统回路模块包括：对象模型模块，鲁棒控制器模块，自适应最优相角陷波滤波器模块；具体步骤为：（1）建立对象模型离散化；（2）设计鲁棒控制器（3）设计最优相角陷波滤波器，并联接入部分系统回路（4）根据P(s)·C(s)的传递函数，采用改进型最小均方值自适应算法，引入H∞混合灵敏度控制设计鲁棒控制器，使陷波中心频率跟踪未知单频干扰频率，从而实现了对随机干扰信号和未知单频干扰信号的最优抑制。可以同时抑制具有宽频随机干扰信号和多个频率未知单频干扰信号等多源干扰信号。

Description

自适应最优相角陷波滤波器设计方法

技术领域

本发明属于控制工程与科学的技术领域，具体涉及一种自适应最优相角陷波滤波器设计方法。

背景技术

在实际所用到的伺服机构中，往往广泛存在着一些集中在窄频带的干扰信号。国内外很多专家学者也对此提出了很多的滤波或者控制方法，比如在低频段(100–600Hz)，常常采用峰值滤波器来有效地抑制低频窄带干扰，然而，峰值滤波器却不能抑制在开环交叉频率周围的中频干扰，因为它的固有相位会损耗相位裕度，在扰动频率附近会使得系统回路灵敏度增益失真。而在中频段(1.6kHz左右)，现有技术常常采用一个相位超前峰值滤波器来抑制系统中出现的中频段窄带干扰，这种滤波器为了保证相位裕度，增加了一个额外的π/2的相位超前的改进微分器将灵敏度曲线平滑化。对于硬盘伺服系统中的高频(4-10kHz)窄带干扰抑制，多采用锁相伺服控制器来抑制由悬架振动引起的偏差干扰。这些滤波器或控制器设计虽然有很好的抑制窄带干扰的能力，但是它们都只能在有限的频率范围有效地抑制窄带干扰。自适应滤波器的出现大大促进了未知频率干扰领域的发展，但是一般的自适应滤波器方法中存在的主要问题是：算法复杂，滤波参数维数过高，自适应时间太长，系统实时性大大降低。

发明内容

本发明针对多个频率变化的未知单频干扰信号和伺服系统本身存在的不确定性等随机干扰在内的多源干扰信号对系统模型的影响，为了实现高精度控制，提供一种自适应最优相角陷波滤波器设计方法。

本发明提供的一种自适应最优相角陷波滤波器设计方法，应用系统回路模块包括：对象模型模块，鲁棒控制器模块和自适应最优相角陷波滤波器模块。鲁棒控制器模块和对象模型模块串联在系统回路的输入端与输出端之间，自适应最优相角陷波滤波器模块并联在系统回路的输入端与鲁棒控制器模块之间，同时串联开关K。

信号传递过程：一方面自适应最优相角陷波滤波器模块的输入信号x₁(n)与误差信号e(n)一起进入自适应最优相角陷波滤波器模块中进行自适应滤波，输出自适应最优相角陷波滤波器模块的输出信号y₁(n)，另一方面x₁(n)信号与y₁(n)信号进行叠加得到鲁棒控制器模块的输入信号x₂(n)，x₂(n)信号经过鲁棒控制器模块和对象模型模块完成控制过程，得到对象模型模块的输出信号y₂(n)，y₂(n)信号与单频干扰信号d(n)叠加，得到误差信号e(n)，误差信号e(n)一方面与x₁(n)信号进入自适应最优相角陷波滤波器模块中，另一方面叠加宽频随机干扰信号，被输入信号r(n)减去，进行负反馈，得到自适应最优相角陷波滤波器模块的输入信号x₁(n)，构成整体的闭环系统。

本发明提供的一种自适应最优相角陷波滤波器设计方法，具体实施步骤包括：

步骤1、通过频率响应实验或物理公式推导建立对象模型模块，用传递函数表示为P(s)，离散化后表示为P(z)。

步骤2、鲁棒控制器模块的传递函数表达式C(s)，离散化后表示为C(z)的设计方法如下：

①设原系统支路灵敏度函数 $S_0\left(s\right)=\frac{1}{1+P\left(s\right)\·C\left(s\right)}=\frac{1}{1+P\left(s\right)C\left(s\right)},$

补灵敏度函数T₀(s)＝1-S₀(s)

②性能加权函数W₁(s)，表示成如下形式：

$W_1\left(s\right)=\frac{K_1}{{[s-(-\mathrm{\ϵ})]}^2}$

K₁为高增益，-ε为二重极点，取值范围为0<ε<1；

③鲁棒加权函数W₂(s)

$\mathrm{\&Delta;}=|\frac{M_{\mathrm{ik}}{e}^{{\mathrm{j\&phi;}}_{\mathrm{ik}}}}{M_i{e}^{{\mathrm{j\&phi;}}_k}}-1|\&le;\left|W_2\left({{\mathrm{j\&omega;}}_i}^{\&prime;}\right)\right|,i=1,...,m;k=1,...,n.$

其中Δ代表对象模型幅值相对误差，j为虚数单位，ω_i'为频率，i=1，…，m，m表示实验重复次数，n表示每次实验的采样点数，k表示第k次实验，(M_ik,φ_ik)为频率ω_i下系统响应的幅值和相角，(M_i,φ_k)为该频率点对应的拟合后的P(s)中的幅值和相角。

W₂(s)的传递函数表达式设计为二阶模型

$W_2\left(s\right)=K_2\frac{s^2+2{\mathrm{\&zeta;}}_1{\mathrm{\&omega;}}_{1}s+{\mathrm{\&omega;}}_1^2}{s^2+2{\mathrm{\&zeta;}}_2{\mathrm{\&omega;}}_{2}s+{\mathrm{\&omega;}}_2^2}$

其中K₂为鲁棒权值函数增益，ζ₁,ζ₂分别为W₂(s)零点多项式和极点多项式的阻尼比，ζ₁,ζ₂∈(0,1)，ω₁,ω₂分别为W₂(s)零点多项式和极点多项式的固有频率；

④得到W₁(s)、W₂(s)、S₀(s)和T₀(s)之后，调用矩阵实验室软件Matrix Laboratory中鲁棒控制器工具箱中的混合灵敏度控制器函数，指令为mixsyn，寻找正则的鲁棒控制器，使得W₁(s)、W₂(s)、S₀(s)和T₀(s)满足无穷范数形式的权值条件：

${\left|\right|\left[\begin{array}{c}{W}_1\left(s\right)S_0\left(s\right)\\ W_2\left(s\right)T_0\left(s\right)\end{array}\right]\left|\right|}_{\&infin;}<\mathrm{\&gamma;}$

其中γ是H_∞的性能指标，取值范围是γ＞0，正则的含义是鲁棒控制器传递函数分母中s的次数要高于分子中s的次数；

得到鲁棒控制器连入闭环系统回路中，既能保证系统的鲁棒稳定性，又能达到抑制宽频域的随机干扰信号的效果。

步骤3、自适应最优相角陷波滤波器模块的传递函数表达式F(s)，离散化后表示为F(z)的设计方法如下：

其中ω₀为需要抑制干扰的频率，s＝jω₀时相角ζ是阻尼比ζ∈(0,1)，K_f是最优相角滤波器正增益。

步骤4、得到z变换之后对象模型模块的传递函数表达式P(z)，鲁棒控制器模块的传递函数表达式C(z)和自适应最优相角陷波滤波器模块的传递函数表达式F(z)以后，系统回路的信号过程如下：

（1）自适应最优相角陷波滤波器模块的输入信号x₁(n)一方面与误差信号e(n)一起进入自适应最优相角陷波滤波器模块中，进行自适应滤波，直到输出自适应最优相角陷波滤波器模块的输出信号y₁(n)；具体步骤如下：

①对于未知频率的单频干扰信号，需要采用根据对象模型改进的最小均方值（LMS）算法，进行部分滤波器参数的自适应过程：

$C\left(z\right)\&CenterDot;P\left(z\right)=\frac{c_0{z}^0+c_1{z}^{-1}+\&CenterDot;\&CenterDot;\&CenterDot;+c_M{z}^{-M}}{1+d_1{z}^{-1}+d_2{z}^{-2}+\&CenterDot;\&CenterDot;\&CenterDot;d_N{z}^{-N}}(M>N)$

c₀,c₁,......,c_M为C(z)·P(z)的分子系数，构成定常矩阵C的元素；M为自然数，d₁,......,d_N为C(z)·P(z)的分母系数，构成定常矩阵D的元素；N为正整数。

②F(s)离散化之后，自适应最优相角陷波滤波器离散化传递函数F(z)在n时刻的表达式为：

$F\left(z\right)=\frac{b_0+b_1{z}^{-1}+b_2{z}^{-2}}{1+a_1\left(n\right)z^{-1}+a_2{z}^{-2}}$

其中b₀，b₁，b₂为F(z)的分子系数，构成定常矩阵B的元素；a₁(n)，a₂为F(z)的分母系数，构成矩阵A(n)的元素；其中A(n)与n时刻有关，a₁(n)随n时刻变化而变化。

③进行自适应滤波后，此时从z⁰到z^-M算起，共M+1个系数分别构成了C(z)·P(z)和F(z)的分子和分母系数矩阵，在n时刻表示为如下形式：

$\left\{\begin{array}{c}A\left(n\right)={[a_1\left(n\right)a_{2}0\&CenterDot;\&CenterDot;\&CenterDot;0]}_{1\&times;(M+1)}\\ B={[b_0{b}_1{b}_{2}0\&CenterDot;\&CenterDot;\&CenterDot;0]}_{1\&times;(M+1)}\\ C={[c_0\&CenterDot;\&CenterDot;\&CenterDot;c_M]}_{1\&times;(M+1)}\\ D={[d_1\&CenterDot;\&CenterDot;\&CenterDot;d_{N}0\&CenterDot;\&CenterDot;\&CenterDot;0]}_{1\&times;(M+1)}\end{array}\right.$

自适应最优相角陷波滤波器模块的输入信号x₁(n)和鲁棒控制器模块的输入信号x2(n)从n时刻往前计算到n-M时刻的值分别为状态矩阵X₁(n)、X₂(n)，自适应最优相角陷波滤波器模块的输出信号y₁(n)和对象模型模块的输出信号y₂(n)从n-1时刻往前计算到n-M-1时刻的值分别为状态矩阵Y₁(n)、Y₂(n)，有如下关系式：

X_i(n)＝[x_i(n)…x_i(n-M)]_1×(M+1),i＝1,2

Y_i(n)＝[y_i(n-1)…y_i(n-M-1)]_1×(M+1),i＝1,2

④由F(z)的表达式得到n时刻的离散状态形式：

$F\left(z\right)=\frac{b_0+b_1{z}^{-1}+b_2{z}^{-2}}{1+a_1\left(n\right)z^{-1}+a_2{z}^{-2}}=\frac{y_1\left(n\right)}{x_1\left(n\right)}$

进而得到：

(1+a₁(n)z^-1+a₂z^-2)y₁(n)＝(b₀+b₁z^-1+b₂z^-2)x₁(n)

又由z^-kx(n)＝x(n-k)，z^-ky(n)＝y(n-k)得

y₁(n)+a₁(n)y₁(n-1)+a₂y₁(n-2)＝b₀x₁(n)+b₁x₁(n-1)+b₂x₁(n-2)

矩阵B与矩阵X₁(n)的转置相结合后，得到结果为

BX₁ ^T(n)=b₀x₁(n)+b₁x₁(n-1)+b₂x₁(n-2)

矩阵A(n)与矩阵Y₁(n)的转置相结合后，得到结果为

A(n)Y₁ ^T(n)=a₁(n)y₁(n-1)+a₂y₁(n-2)

BX₁ ^T(n)的值与A(n)Y₁ ^T(n)的值进行相减操作得到y₁(n)

y₁(n)＝BX₁ ^T(n)-A(n)Y₁ ^T(n)

其中T代表矩阵的转置。

即

y₁(n)＝b₀x₁(n)+b₁x₁(n-1)+b₂x₁(n-2)-a₁(n)y₁(n-1)-a₂y₁(n-2)

（2）x₁(n)信号另一方面与y₁(n)信号进行叠加得到鲁棒控制器模块的输入信号x₂(n)，n时刻x₁(n),x₂(n),y₁(n)数值关系

x₂(n)＝x₁(n)+y₁(n)

从而状态矩阵X₁(n),X₂(n),Y₁(n)有如下关系：

X₂(n)＝X₁(n)+Y₁(n+1)

（3）x₂(n)信号经过鲁棒控制器模块和对象模型模块完成控制过程得到对象模型模块的输出信号y₂(n)，具体操作步骤同得到自适应最优相角陷波滤波器模块的输出信号y₁(n)的步骤相同。

矩阵C与矩阵X₂(n)的转置相结合后，得到结果为

CX₂ ^T(n)=CX₁ ^T(n)+CY₁ ^T(n+1)

选取c_k为矩阵C中按z降幂排序中第一个不为零的系数，0≤k＜N

得到：

CX₁ ^T(n)=c_kx₁(n-k)+…+c_Mx₁(n-M)

CY₁ ^T(n+1)＝c_ky₁(n-k)+…+c_My₁(n-M)

矩阵D与矩阵Y₂(n)的转置相结合后，得到结果为：

DY₂ ^T(n)＝d₁y₂(n-1)-…-d_Ny₂(n-N-1)

由C(z)·P(z)表达式得到离散状态形式

y₂(n)＝CX₂ ^T(n)-DY₂ ^T(n)

y₂(n)=c_kx₁(n-k)+…+c_Mx₁(n-M)+c_ky₁(n-k)+…+c_My₁(n-M)-d₁y₂(n-1)-…-d_Ny₂(n-N-1)

（4）对象模型模块的输出信号y₂(n)与单频干扰信号d(n)叠加，得到误差信号e(n)为：

e(n)＝d(n)+y₂(n)

（5）误差信号e(n)一方面与自适应最优相角陷波滤波器模块的输入信号x₁(n)进入自适应最优相角陷波滤波器模块中，另一方面叠加宽频随机干扰信号，被输入信号r(n)减去，进行负反馈，构成整体的闭环系统。具体步骤如下：

①取误差信号e(n)的均方值为代价函数J(n)：

J(n)＝E[e²(n)]

②代价函数J(n)对于自适应参数a₁(n-k)的梯度为：

$\&dtri;J\left(n\right)=\frac{\&PartialD;\left\{E\right[e^2\left(n\right)\left]\right\}}{\&PartialD;\left[a_1(n-k)\right]},0\&le;k<N$

表示梯度，a₁(n-k)为自适应参数；N为D矩阵中非零元素的个数。

③按照最小均方值算法，采用瞬时值估计梯度矢量：

$\widehat{\&dtri;}J\left(n\right)=\frac{\&PartialD;\left\{E\right[{\hat{e}}^2\left(n\right)\left]\right\}}{\&PartialD;\left[{\hat{a}}_1(n-k)\right]}=2e\left(n\right)\frac{\&PartialD;\left[\hat{e}\left(n\right)\right]}{\&PartialD;\left[{\hat{a}}_1(n-k)\right]},0\&le;k<N$

^表示瞬时值；

④采用改进型LMS迭代算法，求出信号对自适应参数的梯度：

$\frac{\&PartialD;\left[\hat{e}\left(n\right)\right]}{\&PartialD;\left[{\hat{a}}_1(n-k)\right]}=\frac{\&PartialD;\left[\hat{e}\left(n\right)\right]}{\&PartialD;\left[{\hat{y}}_1(n-k)\right]}\&CenterDot;\frac{\&PartialD;\left[{\hat{y}}_1(n-k)\right]}{\&PartialD;\left[{\hat{a}}_1(n-k)\right]}$

⑤由偏导数性质，先求出的偏导数，得到运算结果：

$\frac{\&PartialD;\left[\hat{e}\left(n\right)\right]}{\&PartialD;[{\hat{y}}_1(n-k)}=\frac{\&PartialD;\left[{{\mathrm{CY}}_1}^T(n+1)\right]}{\&PartialD;\left[{\hat{y}}_1(n-k)\right]}=c_k$

⑥由偏导数性质，再求出的偏导数，得到运算结果：

$\frac{\&PartialD;\left[{\hat{y}}_1(n-k)\right]}{\&PartialD;\left[{\hat{a}}_1(n-k)\right]}=-{\hat{y}}_1(n-k-1),0\&le;k<N$

⑦最小均方值算法的瞬时值：

$\widehat{\&dtri;}J\left(n\right)=2\hat{e}\left(n\right)\frac{\&PartialD;\left[\widehat{e\left(n\right)}\right]}{\&PartialD;\left[{\hat{a}}_1(n-k)\right]}=-{2c}_k\hat{e}\left(n\right){\hat{y}}_1(n-k-1),0\&le;k<N$

⑧则自适应参数迭代公式为：

$\widehat{a_1}(n+1)=\widehat{a_1}(n-k)+\frac{(k+1)}{2}\mathrm{\&mu;}[-\widehat{\&dtri;}J\left(n\right)],0\&le;k<N$

μ为迭代收敛步长，满足λ_max是自适应最优相角陷波滤波器模块的输出信号y₁(n)构成的自相关矩阵R的最大特征值。

⑨将的表达式代入上式得到自适应参数迭代公式的最终表达形式：

$\widehat{a_1}(n+1)=\widehat{a_1}(n-k)+(k+1){\mathrm{\&mu;c}}_k\hat{e}\left(n\right){\hat{y}}_1(n-k-1),0\&le;k<N$

将n+1时刻自适应最优相角陷波滤波器系数a₁(n+1)取的值继续运算，到某一个时间ΔT之后，ΔT的大小与收敛步长μ有关，μ越小ΔT越大，达到抑制未知单频干扰信号的效果，自适应最优相角陷波滤波器设计完成。

本发明的优点和积极效果在于：

1)本发明提出了一种自适应最优相角陷波滤波器的设计方法，其滤波器的零点能根据补灵敏度函数在干扰频率处的相角变化而改变，能够在干扰频率处形成陷波，消除已知频率和未知频率的单频干扰信号，设计方法简单有效。

2)本发明研究了一种控制器和滤波器的并联实现结构，进而引入H_∞混合灵敏度控制设计鲁棒控制器，从而实现了对宽频随机干扰信号的最优抑制；

3)本发明基于特殊的最优相角陷波滤波器的模型及其与系统的并联方式，创造性地将其与最小均方值自适应滤波算法结合起来，使陷波滤波器的陷波中心频率能跟踪未知单频干扰信号频率，达到对未知单频干扰信号的抑制作用；

4)本发明在设计最优相角自适应陷波滤波器的过程中，巧妙发现模型中的特点，简化自适应参数为一个，大大缩短自适应时间，增强了实时性，并在一定精度范围内满足抗干扰条件。

5）本发明提出的一种自适应最优相角陷波滤波器设计方法，它能最低限度地影响闭环系统的稳定性，而且在无限频率范围内通过分配零阶滤波器来有效地抑制未知单频干扰信号，并配合设计H_∞混合灵敏度鲁棒控制器，以提高系统的鲁棒性及随机干扰抑制性能，使系统兼具随机干扰抑制能力与未知频率单频干扰信号抑制能力。

附图说明

图1为本发明自适应最优相角陷波滤波器设计方法的流程图。

图2为本发明的自适应最优相角陷波滤波器的数字实现流程图。

具体实施方式

下面将结合附图和具体实施例对本发明作进一步地详细说明。

本发明提出的一种自适应最优相角陷波滤波器设计方法，可以应用于具有宽频随机干扰信号和多个未知频率单频干扰信号等多干扰源信号同时存在的情况下，达到最优抑制。

所述一种自适应最优相角陷波滤波器设计方法所应用的系统回路，如图1所示，所应用的系统回路模块包括：对象模型模块，鲁棒控制器模块以及自适应最优相角陷波滤波器模块。鲁棒控制器模块和对象模型模块串联在系统回路的输入端与输出端之间。自适应最优相角陷波滤波器模块并联在输入端与鲁棒控制器模块之间，同时串联开关K。

鲁棒控制器模块的功能是控制闭环系统鲁棒稳定并且抑制伺服系统本身存在的不确定性相关的随机干扰信号，对象模型模块的功能是用传递函数的方法表述被控对象，自适应最优相角陷波滤波器模块的功能是对未知单频干扰信号的抑制。

首先，在对象数学模型传递函数P(s)基础之上，通过调整性能加权函数W₁(s)和鲁棒加权函数W₂(s)，得出合适的鲁棒控制器传递函数C(s)；其次，设计自适应最优相角陷波滤波器模块的传递函数F(s)，自适应最优相角陷波滤波器模块并联接入部分系统回路中；最后，根据P(s)·C(s)的传递函数，针对滤波器参数a₁采用改进型最小均方值算法，使陷波中心频率跟踪未知单频干扰信号频率。

信号传递方面如图1所示；一方面自适应最优相角陷波滤波器模块的输入信号x₁(n)与误差信号e(n)一起进入自适应最优相角陷波滤波器模块中进行自适应滤波，输出自适应最优相角陷波滤波器模块的输出信号y₁(n)，另一方面x₁(n)信号与y₁(n)信号进行叠加得到鲁棒控制器模块的输入信号x₂(n)，x₂(n)信号经过鲁棒控制器模块和对象模型模块完成控制过程，得到对象模型模块的输出信号y₂(n)，y₂(n)信号与单频干扰信号d(n)叠加，得到误差信号e(n)，误差信号e(n)一方面与x₁(n)信号进入自适应最优相角陷波滤波器模块中，另一方面叠加宽频随机干扰信号，被输入信号r(n)减去，进行负反馈，构成整体的闭环系统。

本发明提供的一种自适应最优相角陷波滤波器设计方法，具体实施步骤包括：

步骤1、通过频率响应实验或物理公式推导建立对象的数学模型，也就是对象模型模块，拉普拉斯变换后对象模型模块的传递函数表达式为P(s)，传递函数P(s)z变换之后得到离散传递函数表达式为P(z)；其中s为拉普拉斯变换算子。

鲁棒控制器模块的传递函数表达式为C(s)，自适应最优相角陷波滤波器模块的传递函数表达式为F(s)，z变换之后鲁棒控制器模块的传递函数表达式为C(z)，自适应最优相角陷波滤波器模块的传递函数表达式为F(z)，z代表z变换的算子。

步骤2、鲁棒控制器模块拉普拉斯变换后的传递函数表达式C(s)，z变换离散化后表示为C(z)的设计方法如下：

①在得到对象模型模块之后，运用相关现代控制理论如H₂和H_∞最优控制技术，使控制对象鲁棒稳定，以H_∞混合灵敏度鲁棒控制器设计为例，对于对象和控制器回路，设原系统支路灵敏度函数为：

$S_0\left(s\right)=\frac{1}{1+P\left(s\right)\&CenterDot;C\left(s\right)}=\frac{1}{1+P\left(s\right)C\left(s\right)}---\left(1\right)$

补灵敏度函数为：

$T_0\left(s\right)=1-S_0\left(s\right)=\frac{P\left(s\right)C\left(s\right)}{1+P\left(s\right)C\left(s\right)}---\left(2\right)$

②W₁(s)为性能加权函数，反映了原系统灵敏度函数在低频段的性能特征及形状要求，从而在低频段获得高增益，增强抗干扰的抑制能力，W₁(s)选择如下形式：

$W_1\left(s\right)=\frac{K_1}{{[s-(-\mathrm{\&epsiv;})]}^2}---\left(3\right)$

K₁为高增益，-ε代表二重极点，取0<ε<1；

③W₂(s)为鲁棒加权函数，因为鲁棒加权函数是不确定函数，由高频段动态和建模参数不确定性所决定，反映了平台系统本身的性能。

如果对象模型模块是稳定的，并且它的传递函数是通过频率响应的实验得到的，测定在一系列的频率下的对象响应的幅值和相角，这一系列频率为ω_i'，i=1，…，m，m表示实验重复次数，另外，实验重复n次，得到的对象模型传递函数P(s)是扫频之后图像拟合而成，扫频的每一个频率点对应的幅值和相角设为(M_ik,φ_ik)，同样输入情况下该频率点对应的拟合后的P(s)图像中的幅值和相角设为(M_i,φ_k)，W₂(s)设计的条件满足：

$\mathrm{\&Delta;}=|\frac{M_{\mathrm{ik}}{e}^{{\mathrm{j\&phi;}}_{\mathrm{ik}}}}{M_i{e}^{{\mathrm{j\&phi;}}_k}}-1|\&le;\left|W_2\left({{\mathrm{j\&omega;}}_i}^{\&prime;}\right)\right|,i=1,...,m;k=1,...,n.---\left(4\right)$

j代表虚数单位，Δ代表对象模型幅值相对误差，n表示每次实验的采样点数，k表示第k次实验。

通过这个方法设计出来的W₂(s)能让|W₂(jω_i')|在所需频率段内将Δ表达式的图线完全包络，其中W₂(s)的传递函数表达式设计为二阶模型：

$W_2\left(s\right)=K_2\frac{s^2+2{\mathrm{\&zeta;}}_1{\mathrm{\&omega;}}_{1}s+{\mathrm{\&omega;}}_1^2}{s^2+2{\mathrm{\&zeta;}}_2{\mathrm{\&omega;}}_{2}s+{\mathrm{\&omega;}}_2^2}---\left(5\right)$

K₂为鲁棒权值函数增益，ζ₁,ζ₂分别为W₂(s)零点多项式和极点多项式的阻尼比，ζ₁,ζ₂∈(0,1)，ω₁,ω₂分别为W₂(s)零点多项式和极点多项式的固有频率，ζ₁ω₁s,ζ₂ω₂s分别代表三个系数乘积；

④得到W₁(s)、W₂(s)、S₀(s)和T₀(s)之后，调用Matrix Laboratory(矩阵实验室软件)中鲁棒控制器工具箱中的混合灵敏度控制器函数，指令为mixsyn，混合灵敏度问题是寻找正则的鲁棒控制器，使得W₁(s)、W₂(s)、S₀(s)和T₀(s)足无穷范数形式的权值条件：

${\left|\right|\left[\begin{array}{c}{W}_1\left(s\right)S_0\left(s\right)\\ W_2\left(s\right)T_0\left(s\right)\end{array}\right]\left|\right|}_{\&infin;}<\mathrm{\&gamma;}---\left(6\right)$

其中γ是H_∞的性能指标，取值范围是γ＞0，正则的含义是鲁棒控制器传递函数分母中s的次数要高于分子中s的次数；

⑤鲁棒控制器模块进行拉普拉斯变换后的传递函数表达式C(s)，z变换离散化之后为C(z)，这样在对数幅频特性曲线中，设计出的原系统支路灵敏度函数S₀(s)在低于交叉频率的频段内幅值小于0dB（分贝），就能达到抑制宽频域的随机干扰信号的效果；

步骤3、自适应最优相角陷波滤波器模块拉普拉斯变换后的传递函数表达式F(s)，z变换离散化后表示为F(z)的设计方法如下：

（1）自适应最优相角陷波滤波器模块的传递函数表达式F(s)，离散化后表示为F(z)的表达式：

其中ω₀为需要抑制干扰的频率，ζ是阻尼比ζ∈(0,1)，K_f是最优相角滤波器正增益，当s＝jω₀时相角

（2）自适应最优相角陷波滤波器模块抑制已知单频干扰信号的过程如下：

①自适应最优相角陷波滤波器模块通过并联结构连入系统中，得出系统总体灵敏度函数S(s)表达式：

$S\left(s\right)=\frac{1}{1+P\left(s\right)C\left(s\right)(1+F\left(s\right))}$ 分子分母同时乘上（1+P(s)C(s)）可得：

$\begin{array}{c}S\left(s\right)=\frac{1+P\left(s\right)C\left(s\right)}{[1+P\left(s\right)C\left(s\right)+P\left(s\right)C\left(s\right)F\left(s\right)]\&CenterDot;(1+P\left(s\right)C\left(s\right))}\\ =\frac{1}{1+P\left(s\right)C\left(s\right)}\&CenterDot;\frac{1+P\left(s\right)C\left(s\right)}{1+P\left(s\right)C\left(s\right)+P\left(s\right)C\left(s\right)F\left(s\right)}\end{array}---\left(8\right)$

打开开关K，从输入到输出的传递函数表示为原系统支路灵敏度函数S₀(s)

$S_0\left(s\right)=\frac{1}{1+P\left(s\right)C\left(s\right)}---\left(9\right)$

合上开关K，滤波器支路从输入到输出的传递函数表示为滤波器支路灵敏度函数S_F(s)，具体地，

最后系统总体灵敏度函数S(s)串行解耦成原系统支路灵敏度函数S₀(s)与滤波器支路灵敏度函数S_F(s)的乘积：

S(s)＝S₀(s)·S_F(s)    （11）

将干扰频率s＝jω₀代入到滤波器支路灵敏度函数S_F(s)中得到相关增益比较关系：

$\left|S_F\left({\mathrm{j\&omega;}}_0\right)\right|=\frac{1}{|1+T_0\left({\mathrm{j\&omega;}}_0\right)F\left({\mathrm{j\&omega;}}_0\right)|}\&GreaterEqual;\frac{1}{1+\left|T_0\left({\mathrm{j\&omega;}}_0\right)\right|\left|F\left({\mathrm{j\&omega;}}_0\right)\right|}---\left(12\right)$

|S_F(jω₀)|的值由|T₀(jω₀)|和|F(jω₀)|决定，由于T₀(jω₀)的值与鲁棒控制器模块C(s)和对象模型模块P(s)有关，已确定，所以只需考虑F(jω₀)的值。

②自适应最优相角陷波滤波器模块的传递函数表达式F(s)，离散化后表示为F(z)的表达式：

其中K_f是最优相角陷波滤波器的正增益，ω₀为需要抑制干扰的频率，ζ是阻尼比ζ∈(0,1)，与补灵敏度函数T₀(s)有关，s＝jω₀时相角两个实零点分别为z₁＝0，

当相角满足如下关系时，滤波器支路灵敏度函数的幅值|S_F(jω₀)|达到最小值：

$\begin{array}{c}F\left({\mathrm{j\&omega;}}_0\right)=K_f\frac{{\mathrm{j\&omega;}}_0\{{\mathrm{\&omega;}}_0\cos \left(\arg \right[T_0\left({\mathrm{j\&omega;}}_0\right)\left]\right)-{\mathrm{j\&omega;}}_0\sin \left(\arg \right[T_0\left({\mathrm{j\&omega;}}_0\right)\left]\right)\}}{{\left({\mathrm{j\&omega;}}_0\right)}^2+2\mathrm{\&zeta;}{\mathrm{\&omega;}}_0\left({\mathrm{j\&omega;}}_0\right)+{{\mathrm{\&omega;}}_0}^2}\\ =K_f\frac{{{\mathrm{j\&omega;}}_0}^2\cos \left(\arg \right[T_0\left({\mathrm{j\&omega;}}_0\right)\left]\right)+{{\mathrm{\&omega;}}_0}^2\sin \left(\arg \right[T_0\left({\mathrm{j\&omega;}}_0\right)\left]\right)}{2\mathrm{j\&zeta;}{{\mathrm{\&omega;}}_0}^2}\\ =K_f\frac{{{\mathrm{\&omega;}}_0}^2\cos \left(\arg \right[T_0\left({\mathrm{j\&omega;}}_0\right)\left]\right)-{{\mathrm{j\&omega;}}_0}^2\sin \left(\arg \right[T_0\left({\mathrm{j\&omega;}}_0\right)\left]\right)}{2\mathrm{\&zeta;}{{\mathrm{\&omega;}}_0}^2}\end{array}---\left(13\right)$

由 $\arg \left[F\left({\mathrm{j\&omega;}}_0\right)\right]=\arctan \frac{\mathrm{Re}\left(F\left({\mathrm{j\&omega;}}_0\right)\right)}{\mathrm{Im}\left(F\left({\mathrm{j\&omega;}}_0\right)\right)}$

Re(F(jω₀))表示对函数F(jω₀)取复数的实部，Im(F(jω₀))表示对函数F(jω₀)取复数的虚部，

可知相角：

$\begin{array}{c}\arg \left[F\left({\mathrm{j\&omega;}}_0\right)\right]=\arctan \frac{\cos \left(\arg \right[T_0\left({\mathrm{j\&omega;}}_0\right)\left]\right)}{-\sin \left(\arg \right[T_0\left({\mathrm{j\&omega;}}_0\right)\left]\right)}=\arctan (-\cot \arg [T_0\left({\mathrm{j\&omega;}}_0\right)\left]\right)=-\arg \left[T_0\left({\mathrm{j\&omega;}}_0\right)\right]\\ \arg \left[T_0\left({\mathrm{j\&omega;}}_0\right)\right]+\arg \left[F\left({\mathrm{j\&omega;}}_0\right)\right]=0\end{array}---\left(14\right)$

所以

$\begin{array}{c}\left|F\left({\mathrm{j\&omega;}}_0\right)\right|=K_f\frac{\sqrt{{\left({{\mathrm{\&omega;}}_0}^2\cos \left(\arg \right[T_0\left({\mathrm{j\&omega;}}_0\right)\left]\right)\right)}^2+{\left({{\mathrm{\&omega;}}_0}^2\sin \left(\arg \right[T_0\left({\mathrm{j\&omega;}}_0\right)\left]\right)\right)}^2}}{\left|2\mathrm{j\&zeta;}{{\mathrm{\&omega;}}_0}^2\right|}\\ =K_f\frac{\sqrt{{{\mathrm{\&omega;}}_0}^4({\cos }^2\left(\arg \right[T_0\left({\mathrm{j\&omega;}}_0\right)\left]\right)+{\sin }^2\left(\arg \right[T_0\left({\mathrm{j\&omega;}}_0\right)\left]\right))}}{\left|2\mathrm{j\&zeta;}{{\mathrm{\&omega;}}_0}^2\right|}\\ =K_f\frac{\sqrt{{{\mathrm{\&omega;}}_0}^4}}{\left|2\mathrm{j\&zeta;}{{\mathrm{\&omega;}}_0}^2\right|}=\frac{K_f}{2\mathrm{\&zeta;}}\end{array}---\left(15\right)$

$|S_F\left({\mathrm{j\&omega;}}_0\right){|}_{\min }=\frac{1}{1+\left|F\left({\mathrm{j\&omega;}}_0\right)\right|\left|T_0\left({\mathrm{j\&omega;}}_0\right)\right|}=\frac{1}{1+\frac{K_f}{2\mathrm{\&zeta;}}\left|T_0\left({\mathrm{j\&omega;}}_0\right)\right|}---\left(16\right)$

K_f与ζ的取值在s＝jω₀的相角满足arg[T₀(jω₀)]+arg[F(jω₀)]＝0时，|S_F(jω₀)|达到最小值，|S_F(jω₀)|_min在对数幅频特性图像上远远低于0dB，进而滤波器支路灵敏度函数S_F(s)与原系统支路灵敏度函数S₀(s)相乘之后，就能在单频干扰信号频率处对原系统支路灵敏度函数S₀(s)产生一个比较大的衰减，并且在其他频率处几乎不改变原系统支路灵敏度函数S₀(s)的形状，而得到最后的系统总体灵敏度函数S(s)，达到抑制已知单频干扰信号的效果。

对于已知频率的单频干扰信号，设计相应的最优相角陷波滤波器来实现抑制作用，而对于未知频率的单频干扰信号，则需要采用根据对象模型模块改进的最小均方值（LMS）算法，对最优相角滤波器进行部分滤波器参数的自适应滤波过程，使得陷波中心频率能跟踪未知单频干扰信号频率，达到对未知频率单频干扰信号的抑制作用。

步骤4、得到z变换之后对象模型模块的传递函数表达式P(z)，鲁棒控制器模块的传递函数表达式C(z)和自适应最优相角陷波滤波器模块的传递函数表达式F(z)以后，系统回路的信号过程如下：

（1）如图1所示：自适应最优相角陷波滤波器模块的输入信号x₁(n)一方面与误差信号e(n)一起进入自适应最优相角陷波滤波器模块中，进行自适应滤波，直到输出自适应最优相角陷波滤波器模块的输出信号y₁(n)，自适应最优相角陷波滤波器的数字实现流程步骤如下：

①改进型最小均方值算法在鲁棒控制器模块及对象模型模块中的实现。为简化算法，考虑鲁棒控制器模块离散化传递函数和对象模型模块离散化传递函数之积：

$C\left(z\right)\&CenterDot;P\left(z\right)=\frac{c_0{z}^0+c_1{z}^{-1}+\&CenterDot;\&CenterDot;\&CenterDot;+c_M{z}^{-M}}{1+d_1{z}^{-1}+d_2{z}^{-2}+\&CenterDot;\&CenterDot;\&CenterDot;d_N{z}^{-N}}(M>N)---\left(17\right)$

其中c₀,c₁,......,c_M为鲁棒控制器离散化传递函数和对象模型离散化传递函数之积C(z)·P(z)的分子系数，该分子系数构成定常矩阵C的元素；M为自然数，d₁,......,d_N为C(z)·P(z)的分母系数，该分母系数构成定常矩阵D的元素；N为正整数。z⁰,......z^-M共M+1个系数为矩阵C，D的列。

②确定自适应参数个数。将自适应最优相角陷波滤波器传递函数离散化之后得到：

$F\left(z\right)=\frac{b_0+b_1{z}^{-1}+b_2{z}^{-2}}{1+a_1{z}^{-1}+a_2{z}^{-2}}---\left(18\right)$

其中，b₀，b₁，b₂代表离散化传递函数F(z)的分子系数，a₁，a₂代表离散化传递函数F(z)的分母系数。

通过实验方法可以验证自适应最优相角陷波滤波器的一个特点：该滤波器零点的变化虽然对滤波器灵敏度函数的形状有一定的影响，但都是在陷波中心频率两边的一些频率段有相应灵敏度增益的增加或者衰减，对于陷波中心频率改变不大，对陷波中心频率处的灵敏度增益也几乎没有改变，即对最优相角陷波滤波器抗干扰性能特征影响不大。所以在一定精度允许范围内，自适应参数可以不考虑滤波器零点参数的相关变化，也就是滤波器分子中的系数b₀，b₁，b₂可以在某一频率认定为不变的，而对于滤波器分母系数a₁，a₂来说，对于只考虑干扰频率ω₀改变的情况下，a₂不变化，也就是说，对于自适应最优相角滤波器算法来说，只需要考虑a₁一个参数的变化。

自适应最优相角陷波滤波器传递函数离散化之后，如图2所示，自适应最优相角陷波滤波器离散化传递函数F(z)在n时刻的表达式为：

$F\left(z\right)=\frac{b_0+b_1{z}^{-1}+b_2{z}^{-2}}{1+a_1\left(n\right)z^{-1}+a_2{z}^{-2}}$

b₀，b₁，b₂构成定常矩阵B的元素；a₁(n)，a₂为F(z)的分母系数，构成矩阵A(n)的元素；其中A(n)与n时刻有关，a₁(n)随n时刻变化而变化。a₁(0)，a₁(1)，......,a₁(k)设为初始值；当n-k时刻得到的自适应参数设为a₁(n-k)；z⁰，……，z^-2为矩阵A(n)，B的列。

③进行自适应滤波后，从z⁰到z^-M算起，共M+1个系数分别构成了C(z)·P(z)和F(z)的分子和分母系数矩阵，其中A(n)与n时刻有关，B、C、D均为定常矩阵，在n时刻表示为如下形式：

$\left\{\begin{array}{c}A\left(n\right)={[a_1\left(n\right)a_{2}0\&CenterDot;\&CenterDot;\&CenterDot;0]}_{1\&times;(M+1)}\\ B={[b_0{b}_1{b}_{2}0\&CenterDot;\&CenterDot;\&CenterDot;0]}_{1\&times;(M+1)}\\ C={[c_0\&CenterDot;\&CenterDot;\&CenterDot;c_M]}_{1\&times;(M+1)}\\ D={[d_1\&CenterDot;\&CenterDot;\&CenterDot;d_{N}0\&CenterDot;\&CenterDot;\&CenterDot;0]}_{1\&times;(M+1)}\end{array}\right.---\left(20\right)$

自适应最优相角陷波滤波器模块的输入信号x₁(n)和鲁棒控制器模块的输入信号x₂(n)从n时刻往前计算到n-M时刻的值分别为状态矩阵分别为X₁(n)和X₂(n)，自适应最优相角陷波滤波器模块的输出信号y₁(n)和对象模型模块的输出信号y₂(n)从n-1时刻往前计算到n-M-1时刻的值分别为状态矩阵Y₁(n)和Y₂(n)，有如下关系式：

X_i(n)＝[x_i(n)…x_i(n-M)]_1×(M+1)(i＝1,2)    (21)

Y_i(n)＝[y_i(n-1)…y_i(n-M-1)]_1×(M+1)(i＝1,2)    (22)

④由F(z)的表达式可得n时刻的离散状态形式：

$F\left(z\right)=\frac{b_0+b_1{z}^{-1}+b_2{z}^{-2}}{1+a_1\left(n\right)z^{-1}+a_2{z}^{-2}}=-\frac{y_1\left(n\right)}{x_1\left(n\right)}--\left(18\right)$

即：

(1+a₁(n)z^-1+a₂z^-2)y₁(n)＝(b₀+b₁z^-1+b₂z^-2)x₁(n)

又由z^-kx(n)＝x(n-k)，z^-ky(n)＝y(n-k)

可得：

y₁(n)+a₁(n)y₁(n-1)+a₂y₁(n-2)＝b₀x₁(n)+b₁x₁(n-1)+b₂x₁(n-2)    （24）

如图2所示，矩阵B与矩阵X₁(n)的转置相结合后，得到结果为

BX₁ ^T(n)=b₀x₁(n)+b₁x₁(n-1)+b₂x₁(n-2)    （25）

矩阵A(n)与矩阵Y₁(n)的转置相结合后，得到结果为

A(n)Y₁ ^T(n)=a₁(n)y₁(n-1)+a₂y₁(n-2)    （26）

⑤BX₁ ^T(n)的值与A(n)Y₁ ^T(n)的值进行相减操作得到y₁(n)

y₁(n)＝BX₁ ^T(n)-A(n)Y₁ ^T(n)    (27)

其中T代表矩阵的转置。

即

y₁(n)＝b₀x₁(n)+b₁x₁(n-1)+b₂x₁(n-2)-a₁(n)y₁(n-1)-a₂y₁(n-2)    （28）

（2）自适应最优相角陷波滤波器模块的输入信号x₁(n)另一方面与自适应最优相角陷波滤波器模块的输出信号y₁(n)进行叠加得到鲁棒控制器模块的输入信号x₂(n)，

x₁(n),x₂(n)和y₁(n)数值关系，如图1所示：

x₂(n)＝x₁(n)+y₁(n)    （29）

从而状态矩阵X₁(n),X₂(n),Y₁(n)有如下关系：

X₂(n)＝X₁(n)+Y₁(n+1)

进行矩阵转置后得到：

X₂ ^T(n)=X₁ ^T(n)+Y₁ ^T(n+1)    （30）

（3）x₂(n)信号经过鲁棒控制器模块和对象模型模块完成控制过程得到对象模型模块的输出信号y₂(n)，得到对象模型模块的输出信号y₂(n)的具体操作步骤，同得到自适应最优相角陷波滤波器模块的输出信号y₁(n)的步骤相同。

矩阵C与矩阵X₂(n)的转置相结合后，得到结果为：

CX₂ ^T(n)=CX₁ ^T(n)+CY₁ ^T(n+1)    （31）

选取c_k为矩阵C中按z降幂排序中第一个不为零的系数，0≤k＜N

得到：

CX₁ ^T(n)=c_kx₁(n-k)+…+c_Mx₁(n-M)    （32）

CY₁ ^T(n+1)＝c_ky₁(n-k)+…+c_My₁(n-M)    （33）

所以：

CX₂ ^T(n)=c_kx₁(n-k)+…+c_Mx₁(n-M)+c_ky₁(n-k)+…+c_My₁(n-M)    （34）

矩阵D与矩阵Y₂(n)的转置相结合后，得到结果为：

DY₂ ^T(n)＝d₁y₂(n-1)-…-d_Ny₂(n-N-1)    （35）

由C(z)·P(z)表达式可得到离散状态形式：

y₂(n)=CX₂ ^T(n)-DY₂ ^T(n)    （36）

y₂(n)＝c_kx₁(n-k)+…+c_Mx₁(n-M)+c_ky₁(n-k)+…+c_My₁(n-M)-d₁y₂(n-1)-…-d_Ny₂(n-N-1)（37）

（4）对象模型模块的输出信号y₂(n)与单频干扰信号d(n)叠加，得到误差信号e(n)，如图1所示：n时刻的误差信号e(n)可表示为：

（5）误差信号e(n)一方面与自适应最优相角陷波滤波器模块的输入信号x₁(n)进入自适应最优相角陷波滤波器模块中，进行自适应滤波过程。另一方面叠加宽频随机干扰，被输入信号r(n)减去，进行负反馈，构成整体的闭环系统。自适应滤波过程步骤如下：

①设误差信号e(n)的均方值代价函数为J(n)：

J(n)＝E[e²(n)]    （39）

②代价函数J(n)对于自适应参数a₁(n-k)的梯度为：

$\&dtri;J\left(n\right)=\frac{\&PartialD;\left\{E\right[e^2\left(n\right)\left]\right\}}{\&PartialD;\left[a_1(n-k)\right]},0\&le;k<N---\left(40\right)$

该改进算法的目的即满足代价函数J(n)趋向最小值，采用改进型LMS迭代算法，在常规的梯度迭代算法中一般要考虑延时最少相关的误差信号与迭代参数，即n时刻的误差信号对n-q时刻的迭代参数求偏导，用表示梯度，其中q越小收敛效果越好，设a₁(n-k)为自适

应参数；N为D矩阵中非零元素的个数。

③按照最小均方值算法，采用瞬时值估计梯度矢量的方法：

$\widehat{\&dtri;}J\left(n\right)=\frac{\&PartialD;\left\{E\right[\widehat{e^2\left(n\right)}]}{\&PartialD;\left[{\hat{a}}_1(n-k)\right]}=2e\left(n\right)\frac{\&PartialD;\left[\widehat{e\left(n\right)}\right]}{\&PartialD;\left[{\hat{a}}_1(n-k)\right]}0\&le;k<N---\left(41\right)$

^表示瞬时值；

④信号对自适应参数的梯度计算过程如下：

对误差信号e(n)，CY₁ ^T(n+1)与自适应参数a₁(n-k)相关，又由于CX₁ ^T(n)与a₁(n-k)无关，DY₂ ^T(n)与a₁(n-k-1)有关，即e(n)求偏导时只用考虑CY₁ ^T(n+1)对a₁(n-k)求偏导，由于CY₁ ^T(n+1)中与a₁(n-k)有关的项只有c_ky₁(n-k)，即：c_ky₁(n-k)对y₁(n-k)的梯度计算过程：

由n时刻的e(n)的展开式，与a₁(n-k)有关的项只有y₁(n-k)，故由二阶偏微分性质可得：

$\frac{\&PartialD;\left[\widehat{e\left(n\right)}\right]}{\&PartialD;\left[\widehat{a_1}(n-k)\right]}=\frac{\&PartialD;\left[\widehat{e\left(n\right)}\right]}{\&PartialD;\left[{\hat{y}}_1(n-k)\right]}\&CenterDot;\frac{\&PartialD;[{\hat{y}}_1(n-k)}{\&PartialD;[{\hat{a}}_1(n-k)}---\left(42\right)$

⑤由偏导数性质，先求出对的偏导数，得到运算结果

$\frac{\&PartialD;\left[\widehat{e\left(n\right)}\right]}{\&PartialD;\left[{\hat{y}}_1(n-k)\right]}=\frac{\&PartialD;[d\left(n\right)+c_k{x}_1(n-k)+\&CenterDot;\&CenterDot;\&CenterDot;+c_M{x}_1(n-M)+c_k{y}_1(n-k)+\&CenterDot;\&CenterDot;\&CenterDot;+c_M{y}_1(n-M)-(d_1{y}_2(n-1)-\&CenterDot;\&CenterDot;\&CenterDot;-d_N{y}_2(n-N-1)]}{\&PartialD;\left[{\hat{y}}_1(n-k)\right]}$

$\frac{\&PartialD;\left[\widehat{e\left(n\right)}\right]}{\&PartialD;\left[{\hat{y}}_1(n-k)\right]}=\frac{\&PartialD;\left[{{\mathrm{CY}}_1}^T(n+1)\right]}{\&PartialD;\left[{\hat{y}}_1(n-k)\right]}=c_k---\left(43\right)$

⑥自适应最优相角陷波滤波器模块的输出信号y₁(n)，经过k+1次延时环节之后，得到y₁(n-k-1)，

如附图2所示，结合公式z^-ka₁(n)＝a₁(n-k)可得：

z^-(k+1)y₁(n)＝y₁(n-k-1)

由偏导数性质，再求出对的偏导数，得到运算结果：

$\begin{array}{c}\frac{\&PartialD;\left[{\hat{y}}_1(n-k)\right]}{\&PartialD;\left[{\hat{a}}_1(n-k)\right]}=\frac{\&PartialD;[b_0{\hat{x}}_1(n-k)+b_1{\hat{x}}_1(n-k-1)+b_2{\hat{x}}_1(n-k-2)-\widehat{a_1}(n-k){\hat{y}}_1(n-k-1)-a_2{\hat{y}}_1(n-k-2)]}{\&PartialD;\left[{\hat{a}}_1(n-k)\right]}\\ =-\widehat{y_1}(n-k-1),0\&le;k<N\end{array}---\left(44\right)$

⑦最小均方值算法的瞬时值：

$\widehat{\&dtri;}J\left(n\right)=2\hat{e}\left(n\right)\frac{\&PartialD;\left[\widehat{e\left(n\right)}\right]}{\&PartialD;\left[{\hat{a}}_1(n-k)\right]}=-{2c}_k\hat{e}\left(n\right){\hat{y}}_1(n-k-1),0\&le;k<N---\left(45\right)$

⑧如附图2所示，得到y₁(n)后，设定前k+1个时刻的值a₁(0)，a₁(1)，......,a₁(k)为初始值，从第k+1时刻开始，每隔k+1进行一次迭代计算，当第n-k时刻的自适应参数a₁(n-k)进行迭代计算后，对应的自适应参数值为a₁(n+1)；

自适应参数a₁(n-k)与k+1倍的单次参数修正值相加，得到n+1时刻的自适应参数值a₁(n+1)，其中单次参数修正值为迭代收敛步长μ乘以二分之一倍的负的代价函数的梯度，则自适应参数迭代公式为：

$\widehat{a_1}(n+1)=\widehat{a_1}(n-k)+\frac{(k+1)}{2}\mathrm{\&mu;}[-\frac{1}{2}\widehat{\&dtri;}J\left(n\right)],0\&le;k<N---\left(46\right)$

其中，μ为迭代收敛步长，满足λ_max是输出信号y₁(n)构成的自相关矩阵R的最大特征值。

⑨将的表达式代入上式可得自适应参数迭代公式的最终表达形式：

$\widehat{a_1}(n+1)=\widehat{a_1}(n-k)+(k+1){\mathrm{\&mu;c}}_k\hat{e}\left(n\right){\hat{y}}_1(n-k-1),0\&le;k<N---\left(47\right)$

将n+1时刻自适应最优相角陷波滤波器系数a₁(n+1)取的值继续运算，到某一个时间ΔT之后，ΔT的大小与收敛步长μ有关，μ越小，ΔT越大，自适应最优相角陷波滤波器的陷波中心频率会跟踪上单频干扰信号频率，自适应过程完成，从而达到抑制单频干扰信号的效果，自适应最优相角陷波滤波器设计完成。

该设计方法重点研究了一种特殊的最优相角陷波滤波器的结构和自适应算法，并以此为基础研究了一种控制器和滤波器的并联实现结构，进而引入H_∞混合灵敏度控制设计鲁棒控制器，从而实现了对随机干扰信号和未知单频干扰信号信号的最优抑制。这一算法可以应用于具有宽频随机干扰，多个未知频率的单频干扰信号等多源干扰同时存在情况的最优抑制。

本发明说明书中未作详细描述的内容属于本领域专业技术人员公知的现有技术。

Claims (3)

1.自适应最优相角陷波滤波器设计方法，其特征在于：系统回路模块包括：对象模型模块，鲁棒控制器模块，自适应最优相角陷波滤波器模块；

鲁棒控制器模块和对象模型模块串联在系统回路的输入端与输出端之间；自适应最优相角陷波滤波器模块并联在系统回路的输入端与鲁棒控制器模块之间，同时串联开关K；

自适应最优相角陷波滤波器设计方法的具体实施步骤包括：

步骤1、建立对象模型模块，拉普拉斯变换后传递函数表示为P(s)，z变换离散化后表示为P(z)；s为拉普拉斯变换算子；

步骤2、鲁棒控制器模块拉普拉斯变换后的传递函数表达式C(s)，z变换离散化后表示为C(z)的设计方法如下：

①设原系统支路灵敏度函数 $S_0\left(s\right)=\frac{1}{1+P\left(s\right)\&CenterDot;C\left(s\right)}=\frac{1}{1+P\left(s\right)C\left(s\right)},$

补灵敏度函数T₀(s)＝1-S₀(s)

②性能加权函数W₁(s)，表示成如下形式：

$W_1\left(s\right)=\frac{K_1}{{[s-(-\mathrm{\&epsiv;})]}^2}$

K₁为高增益，-ε为二重极点，取值范围为0<ε<1；

③鲁棒加权函数W₂(s)的传递函数表达式设计为二阶模型

$W_2\left(s\right)=K_2\frac{s^2+2{\mathrm{\&zeta;}}_1{\mathrm{\&omega;}}_{1}s+{\mathrm{\&omega;}}_1^2}{s^2+2{\mathrm{\&zeta;}}_2{\mathrm{\&omega;}}_{2}s+{\mathrm{\&omega;}}_2^2}$

其中K₂为鲁棒权值函数增益，ζ₁,ζ₂分别为W₂(s)零点多项式和极点多项式的阻尼比，ζ₁,ζ₂∈(0,1)，ω₁,ω₂分别为W₂(s)零点多项式和极点多项式的固有频率；

④得到W₁(s)、W₂(s)、S₀(s)和T₀(s)之后，调用矩阵实验室软件Matrix Laboratory中鲁棒控制器工具箱中的混合灵敏度控制器函数，指令为mixsyn，混合灵敏度问题是寻找正则的鲁棒控制器，使得W₁(s)、W₂(s)、S₀(s)和T₀(s)满足无穷范数形式的权值条件：

${\left|\right|\left[\begin{array}{c}{W}_1\left(s\right)S_0\left(s\right)\\ W_2\left(s\right)T_0\left(s\right)\end{array}\right]\left|\right|}_{\&infin;}<\mathrm{\&gamma;}$

其中γ是H_∞的性能指标，取值范围是γ＞0；将得到的正则的鲁棒控制器连入闭环系统回路中；

步骤3、自适应最优相角陷波滤波器模块的传递函数表达式F(s)，离散化后表示为F(z)的设计方法如下：

其中ω₀为需要抑制干扰的频率，s＝jω₀时相角ζ是阻尼比ζ∈(0,1)，K_f是最优相角滤波器正增益；

步骤4、得到z变换之后对象模型模块的传递函数表达式P(z)，鲁棒控制器模块的传递函数表达式C(z)和自适应最优相角陷波滤波器模块的传递函数表达式F(z)以后，系统回路的信号过程如下：

（1）自适应最优相角陷波滤波器模块的输入信号x₁(n)与误差信号e(n)一起进入自适应最优相角陷波滤波器模块中，进行自适应滤波，直到得到自适应最优相角陷波滤波器模块的输出信号y₁(n)；具体步骤如下：

①对于未知频率的单频干扰信号，需要采用根据对象模型改进的最小均方值算法，进行部分滤波器参数的自适应过程：

$C\left(z\right)\&CenterDot;P\left(z\right)=\frac{c_0{z}^0+c_1{z}^{-1}+\&CenterDot;\&CenterDot;\&CenterDot;+c_M{z}^{-M}}{1+d_1{z}^{-1}+d_2{z}^{-2}+\&CenterDot;\&CenterDot;\&CenterDot;d_N{z}^{-N}}(M>N)$

c₀,c₁,......,c_M为C(z)·P(z)的分子系数，构成定常矩阵C的元素；M为自然数，d₁,......,d_N为C(z)·P(z)的分母系数，构成定常矩阵D的元素；N为正整数；

②F(s)离散化之后，自适应最优相角陷波滤波器离散化传递函数F(z)在n时刻的表达式为：

$F\left(z\right)=\frac{b_0+b_1{z}^{-1}+b_2{z}^{-2}}{1+a_1\left(n\right)z^{-1}+a_2{z}^{-2}}$

其中b₀，b₁，b₂为F(z)的分子系数，构成定常矩阵B的元素；a₁(n)，a₂为F(z)的分母系数，构成矩阵A(n)的元素；其中A(n)与n时刻有关，a₁(n)随n时刻变化而变化；

③进行自适应滤波后，此时从z⁰到z^-M算起，共M+1个系数分别构成了C(z)·P(z)和F(z)的分子和分母系数矩阵，在n时刻表示为如下形式：

$\left\{\begin{array}{c}A\left(n\right)={[a_1\left(n\right)a_{2}0\&CenterDot;\&CenterDot;\&CenterDot;0]}_{1\&times;(M+1)}\\ B={[b_0{b}_1{b}_{2}0\&CenterDot;\&CenterDot;\&CenterDot;0]}_{1\&times;(M+1)}\\ C={[c_0\&CenterDot;\&CenterDot;\&CenterDot;c_M]}_{1\&times;(M+1)}\\ D={[d_1\&CenterDot;\&CenterDot;\&CenterDot;d_{N}0\&CenterDot;\&CenterDot;\&CenterDot;0]}_{1\&times;(M+1)}\end{array}\right.$

自适应最优相角陷波滤波器模块的输入信号x₁(n)和鲁棒控制器模块的输入信号x₂(n)从n时刻往前计算到n-M时刻的值分别为状态矩阵X₁(n)、X₂(n)，自适应最优相角陷波滤波器模块的输出信号y₁(n)和对象模型模块的输出信号y₂(n)从n-1时刻往前计算到n-M-1时刻的值分别为状态矩阵Y₁(n)、Y₂(n)，有如下关系式：

X_i(n)＝[x_i(n)…x_i(n-M)]_1×(M+1)i＝1,2

Y_i(n)＝[y_i(n-1)…y_i(n-M-1)]_1×(M+1),i＝1,2

④由F(z)的表达式得到n时刻的离散状态形式：

$F\left(z\right)=\frac{b_0+b_1{z}^{-1}+b_2{z}^{-2}}{1+a_1\left(n\right)z^{-1}+a_2{z}^{-2}}=\frac{y_1\left(n\right)}{x_1\left(n\right)}$

(1+a₁(n)z^-1+a₂z^-2)y₁(n)＝(b₀+b₁z^-1+b₂z^-2)x₁(n)

结合z^-kx(n)＝x(n-k)，z^-ky(n)＝y(n-k)得：

y₁(n)+a₁(n)y₁(n-1)+a₂y₁(n-2)＝b₀x₁(n)+b₁x₁(n-1)+b₂x₁(n-2)

矩阵B与矩阵X₁(n)的转置相结合后，得到结果为

BX₁ ^T(n)=b₀x₁(n)+b₁x₁(n-1)+b₂x₁(n-2)

矩阵A(n)与矩阵Y₁(n)的转置相结合后，得到结果为

A(n)Y₁ ^T(n)=a₁(n)y₁(n-1)+a₂y₁(n-2)

BX₁ ^T(n)的值与A(n)Y₁ ^T(n)的值进行相减操作得到y₁(n)

y₁(n)＝BX₁ ^T(n)-A(n)Y₁ ^T(n)

其中T代表矩阵的转置，

y₁(n)＝b₀x₁(n)+b₁x₁(n-1)+b₂x₁(n-2)-a₁(n)y₁(n-1)-a₂y₁(n-2)

（2）x₁(n)信号另一方面与y₁(n)信号进行叠加得到鲁棒控制器模块的输入信号x₂(n)，n时刻x₁(n),x₂(n),y₁(n)数值关系：

x₂(n)＝x₁(n)+y₁(n)

从而状态矩阵X₁(n),X₂(n),Y₁(n)有如下关系：

X₂(n)＝X₁(n)+Y₁(n+1)

（3）x₂(n)信号经过鲁棒控制器模块和对象模型模块完成控制过程得到对象模型模块的输出信号y₂(n)，具体操作步骤同得到自适应最优相角陷波滤波器模块的输出信号y₁(n)的步骤相同；

矩阵C与矩阵X₂(n)的转置相结合后，得到结果为

CX₂ ^T(n)=CX₁ ^T(n)+CY₁ ^T(n+1)

选取c_k为矩阵C中按z降幂排序中第一个不为零的系数，0≤k＜N

得到：

CX₁ ^T(n)=c_kx₁(n-k)+…+c_Mx₁(n-M)

CY₁ ^T(n+1)＝c_ky₁(n-k)+…+c_My₁(n-M)

矩阵D与矩阵Y₂(n)的转置相结合后，得到结果为：

DY₂ ^T(n)＝d₁y₂(n-1)-…-d_Ny₂(n-N-1)

由C(z)·P(z)表达式得到离散状态形式

y₂(n)＝CX₂ ^T(n)-DY₂ ^T(n)

y₂(n)=c_kx₁(n-k)+…+c_Mx₁(n-M)+c_ky₁(n-k)+…+c_My₁(n-M)-d₁y₂(n-1)-…-d_Ny₂(n-N-1)

（4）对象模型模块的输出信号y₂(n)与单频干扰信号d(n)叠加，得到误差信号e(n)：

e(n)＝d(n)+y₂(n)

（5）误差信号e(n)一方面与自适应最优相角陷波滤波器模块的输入信号x₁(n)进入自适应最优相角陷波滤波器模块中，另一方面叠加宽频随机干扰信号，被输入信号r(n)减去，进行负反馈，构成整体的闭环系统；具体步骤如下：

①取误差信号e(n)的均方值为代价函数J(n)：

J(n)＝E[e²(n)]

②代价函数J(n)对于自适应参数a₁(n-k)的梯度为：

$\&dtri;J\left(n\right)=\frac{\&PartialD;\left\{E\right[e^2\left(n\right)\left]\right\}}{\&PartialD;\left[a_1(n-k)\right]},0\&le;k<N$

表示梯度，a₁(n-k)为自适应参数；

③按照最小均方值算法，采用瞬时值估计梯度矢量：

$\widehat{\&dtri;}J\left(n\right)=\frac{\&PartialD;\left\{E\right[{\hat{e}}^2\left(n\right)\left]\right\}}{\&PartialD;\left[{\hat{a}}_1(n-k)\right]}=2e\left(n\right)\frac{\&PartialD;\left[\hat{e}\left(n\right)\right]}{\&PartialD;\left[{\hat{a}}_1(n-k)\right]},0\&le;k<N$

^表示瞬时值；

④确定信号对自适应参数的梯度：

$\frac{\&PartialD;\left[\hat{e}\left(n\right)\right]}{\&PartialD;\left[{\hat{a}}_1(n-k)\right]}=\frac{\&PartialD;\left[\hat{e}\left(n\right)\right]}{\&PartialD;\left[{\hat{y}}_1(n-k)\right]}\&CenterDot;\frac{\&PartialD;\left[{\hat{y}}_1(n-k)\right]}{\&PartialD;\left[{\hat{a}}_1(n-k)\right]}$

⑤由偏导数性质，先求出对的偏导数，得到：

$\frac{\&PartialD;\left[\hat{e}\left(n\right)\right]}{\&PartialD;[{\hat{y}}_1(n-k)}=\frac{\&PartialD;\left[{{\mathrm{CY}}_1}^T(n+1)\right]}{\&PartialD;\left[{\hat{y}}_1(n-k)\right]}=c_k$

⑥由偏导数性质，再求出对的偏导数，得到：

$\frac{\&PartialD;\left[{\hat{y}}_1(n-k)\right]}{\&PartialD;\left[{\hat{a}}_1(n-k)\right]}=-{\hat{y}}_1(n-k-1),0\&le;k<N$

⑦最小均方值算法的瞬时值：

$\widehat{\&dtri;}J\left(n\right)=2\hat{e}\left(n\right)\frac{\&PartialD;\left[\widehat{e\left(n\right)}\right]}{\&PartialD;\left[{\hat{a}}_1(n-k)\right]}=-{2c}_k\hat{e}\left(n\right){\hat{y}}_1(n-k-1),0\&le;k<N$

⑧则自适应参数迭代公式为：

$\widehat{a_1}(n+1)=\widehat{a_1}(n-k)+\frac{(k+1)}{2}\mathrm{\&mu;}[-\widehat{\&dtri;}J\left(n\right)],0\&le;k<N$

μ为迭代收敛步长，满足λ_max是自适应最优相角陷波滤波器模块的输出信号y₁(n)构成的自相关矩阵的最大特征值；

⑨获得自适应参数迭代公式的最终表达形式：

$\widehat{a_1}(n+1)=\widehat{a_1}(n-k)+(k+1){\mathrm{\&mu;c}}_k\hat{e}\left(n\right){\hat{y}}_1(n-k-1),0\&le;k<N$

将n+1时刻自适应最优相角陷波滤波器系数a₁(n+1)取的值进行运算，到某一个时间ΔT之后，μ越小ΔT越大，达到抑制未知单频干扰信号的效果，自适应最优相角陷波滤波器设计完成。

2.如权利要求1所述的自适应最优相角陷波滤波器设计方法，其特征在于：鲁棒加权函数W₂(s)通过如下方法设计：

$\mathrm{\&Delta;}=|\frac{M_{\mathrm{ik}}{e}^{{\mathrm{j\&phi;}}_{\mathrm{ik}}}}{M_i{e}^{{\mathrm{j\&phi;}}_k}}-1|\&le;\left|W_2\left({{\mathrm{j\&omega;}}_i}^{\&prime;}\right)\right|,i=1,...,m;k=1,...,n.$

其中Δ代表对象模型幅值相对误差，j为虚数单位，ω_i'为频率，i=1，…，m，m表示实验重复次数，n表示每次实验的采样点数，k表示第k次实验，(M_ik,φ_ik)为频率ω_i下系统响应的幅值和相角，(M_i,φ_k)为该频率点对应的拟合后的P(s)中的幅值和相角。

3.如权利要求1所述的自适应最优相角陷波滤波器设计方法，其特征在于：自适应最优相角陷波滤波器模块抑制已知单频干扰信号的过程如下：

①自适应最优相角陷波滤波器模块通过并联结构连入系统中，得出系统总体灵敏度函数S(s)表达式：

$\begin{array}{c}S\left(s\right)=\frac{1}{1+P\left(s\right)C\left(s\right)(1+F\left(s\right))}\\ =\frac{1}{1+P\left(s\right)C\left(s\right)}\&CenterDot;\frac{1+P\left(s\right)C\left(s\right)}{1+P\left(s\right)C\left(s\right)+P\left(s\right)C\left(s\right)F\left(s\right)}\end{array}$

打开开关K，从输入到输出的传递函数表示为原系统支路灵敏度函数S₀(s)：

$S_0\left(s\right)=\frac{1}{1+P\left(s\right)C\left(s\right)}$

合上开关K，滤波器支路从输入到输出的传递函数表示为滤波器支路灵敏度函数S_F(s)：

$\begin{array}{c}{S}_F\left(s\right)=\frac{1+P\left(s\right)C\left(s\right)}{1+P\left(s\right)C\left(s\right)+P\left(s\right)C\left(s\right)F\left(s\right)}\\ =\frac{1}{1+\frac{P\left(s\right)C\left(s\right)}{1+P\left(s\right)C\left(s\right)}F\left(s\right)}\\ =\frac{1}{1+F\left(s\right)T_0\left(s\right)}\end{array}$

最后系统总体灵敏度函数S(s)串行解耦成原系统支路灵敏度函数S₀(s)与滤波器支路灵敏度函数S_F(s)的乘积：

S(s)＝S₀(s)·S_F(s)

将干扰频率s＝jω₀代入到滤波器支路灵敏度函数S_F(s)中得到相关增益比较关系：

$\left|S_F\left({\mathrm{j\&omega;}}_0\right)\right|=\frac{1}{|1+T_0\left({\mathrm{j\&omega;}}_0\right)F\left({\mathrm{j\&omega;}}_0\right)|}\&GreaterEqual;\frac{1}{1+\left|T_0\left({\mathrm{j\&omega;}}_0\right)\right|\left|F\left({\mathrm{j\&omega;}}_0\right)\right|}$

②自适应最优相角陷波滤波器模块的传递函数表达式F(s)，离散化后表示为F(z)；

当相角满足如下关系时，滤波器支路灵敏度函数的幅值|S_F(jω₀)|达到最小值：

$\begin{array}{c}F\left({\mathrm{j\&omega;}}_0\right)=K_f\frac{{\mathrm{j\&omega;}}_0\{{\mathrm{\&omega;}}_0\cos \left(\arg \right[T_0\left({\mathrm{j\&omega;}}_0\right)\left]\right)-{\mathrm{j\&omega;}}_0\sin \left(\arg \right[T_0\left({\mathrm{j\&omega;}}_0\right)\left]\right)\}}{{\left({\mathrm{j\&omega;}}_0\right)}^2+2\mathrm{\&zeta;}{\mathrm{\&omega;}}_0\left({\mathrm{j\&omega;}}_0\right)+{{\mathrm{\&omega;}}_0}^2}\\ =K_f\frac{{{\mathrm{\&omega;}}_0}^2\cos \left(\arg \right[T_0\left({\mathrm{j\&omega;}}_0\right)\left]\right)-{{\mathrm{j\&omega;}}_0}^2\sin \left(\arg \right[T_0\left({\mathrm{j\&omega;}}_0\right)\left]\right)}{2\mathrm{\&zeta;}{{\mathrm{\&omega;}}_0}^2}\end{array}$

相角 $\begin{array}{c}\arg \left[F\left({\mathrm{j\&omega;}}_0\right)\right]=\arctan \frac{\cos \left(\arg \right[T_0\left({\mathrm{j\&omega;}}_0\right)\left]\right)}{-\sin \left(\arg \right[T_0\left({\mathrm{j\&omega;}}_0\right)\left]\right)}=\arctan (-\cot \arg [T_0\left({\mathrm{j\&omega;}}_0\right)\left]\right)=-\arg \left[T_0\left({\mathrm{j\&omega;}}_0\right)\right]\\ \arg \left[T_0\left({\mathrm{j\&omega;}}_0\right)\right]+\arg \left[F\left({\mathrm{j\&omega;}}_0\right)\right]=0\end{array}$

所以

$\begin{array}{c}\left|F\left({\mathrm{j\&omega;}}_0\right)\right|=K_f\frac{\sqrt{{\left({{\mathrm{\&omega;}}_0}^2\cos \left(\arg \right[T_0\left({\mathrm{j\&omega;}}_0\right)\left]\right)\right)}^2+{\left({{\mathrm{\&omega;}}_0}^2\sin \left(\arg \right[T_0\left({\mathrm{j\&omega;}}_0\right)\left]\right)\right)}^2}}{\left|2\mathrm{j\&zeta;}{{\mathrm{\&omega;}}_0}^2\right|}\\ =K_f\frac{\sqrt{{{\mathrm{\&omega;}}_0}^4}}{|2\mathrm{j\&zeta;}{{\mathrm{\&omega;}}_0}^2}=\frac{K_f}{2\mathrm{\&zeta;}}\end{array}$

$|S_F\left({\mathrm{j\&omega;}}_0\right){|}_{\min }=\frac{1}{1+\left|F\left({\mathrm{j\&omega;}}_0\right)\right|\left|T_0\left({\mathrm{j\&omega;}}_0\right)\right|}=\frac{1}{1+\frac{K_f}{2\mathrm{\&zeta;}}\left|T_0\left({\mathrm{j\&omega;}}_0\right)\right|}$

K_f与ζ的取值在s＝jω₀的相角满足arg[T₀(jω₀)]+arg[F(jω₀)]＝0时，|S_F(jω₀)|达到最小值，进而滤波器支路灵敏度函数S_F(s)与原系统支路灵敏度函数S₀(s)相乘之后，在单频干扰信号频率处对原系统支路灵敏度函数S₀(s)产生比较大的衰减，并且在其他频率处几乎不改变原系统支路灵敏度函数S₀(s)的形状，得到系统总体灵敏度函数S(s)，达到抑制已知单频干扰信号的效果。