CN105783898A - 一种基于频域自适应lms算法的磁悬浮转子谐波振动抑制方法 - Google Patents

Abstract

本发明公开了一种基于频域自适应LMS算法的磁悬浮转子谐波振动抑制方法。本发明主要针对四自由度主被动磁悬浮控制力矩陀螺，在考虑转子不平衡质量和传感器谐波的情况下，对磁悬浮转子进行动力学建模，然后设计了一种基于频域自适应LMS(Least Mean Square)算法的谐波振动抑制方法对转子系统的倍频振动进行抑制，并使用陷波器进行同频振动的抑制。本发明能够抑制磁轴承线圈电流中的谐波分量，进而抑制转子位移中的谐波分量，适用于存在质量不平衡和传感器谐波的磁悬浮转子系统谐波振动的主动控制。

Description

一种基于频域自适应LMS算法的磁悬浮转子谐波振动抑制方法

技术领域

本发明是一种基于频域自适应LMS(LeastMeanSquare)算法的磁悬浮转子谐波振动力抑制方法，针对磁悬浮控制力矩陀螺的谐波振动进行抑制，从而使以磁悬浮控制力矩陀螺为惯性执行机构的卫星平台达到“超稳超静”的要求。

背景技术

高分辨率对地观测，卫星间激光通信，都需要高指向精度、高姿态稳定度的卫星平台提供可靠、微振动的“超静”条件，所以，对超静卫星平台的研究具有非常重要的现实意义。相比传统的机械轴承惯性执行机构，磁悬浮控制力矩陀螺采用磁轴承支承，无摩擦，长寿命，并且能够进行主动控制，从而实现微振动，是实现“超静”卫星平台的理想执行机构之一。

根据磁悬浮转子可控自由度的数目，可将磁悬浮CMG(ControlMomentGyro)分为全主动磁悬浮CMG和主被动磁悬浮CMG两类。全主动磁悬浮CMG除电机驱动自由度之外，其余的五个自由度全部由主动磁轴承来实现稳定的悬浮。而主被动磁悬浮CMG除电机驱动自由度之外，其余的五个自由度不完全主动可控，部分自由度由被动磁轴承实现被动稳定悬浮，无需主动控制的参与。虽然理论上磁悬浮CMG无摩擦和振动，但是由于转子材料密度不均匀，加工误差，传感器装配误差以及电磁不均匀等因素的影响，仍然存在高频振动源，转子系统受到多谐波振动力的影响。振动源主要有转子自身的质量不平衡和位移传感器的测量噪声，即传感器谐波。当谐波频率接近或达到转子框架或CMG壳体的模态时，会引发共振，导致系统失稳，因此需要对谐波进行抑制。

现有的谐波抑制的方法主要分为两类，一类是通过串联不同频率的陷波器分别对各个谐波进行抑制；另一类通过估计转子不平衡和传感器谐波，采用自适应控制使设定的目标函数收敛。多个陷波器串联的方法直接导致算法的计算量很大，而且需要考虑不同滤波器间的收敛问题，设计较为复杂，难以满足实时性、快速性等要求。相比传统的LMS算法，频域LMS算法计算量更小，通过频域完成权值向量的自适应更新，并且无需串联多个滤波器，一个滤波器可以同时对不同频率的扰动进行抑制。现有技术存在以下问题：(1)目前的谐波振动抑制算法计算量较大，实时性受到限制。(2)使用频域LMS算法进行谐波抑制需要权衡算法的收敛速度和稳态精度。

发明内容

本发明的目的：针对目前已有技术的不足，发明一种基于频域自适应LMS算法的磁悬浮转子谐波振动控制方法。利用快速傅里叶变换算法在频域上完成滤波器系数的自适应，采用50％的重叠存储法使运算效率达到最高，有效降低了算法的计算量，并且通过实时的改变步长和块长提高算法的性能。

本发明的技术解决方案：一种基于频域自适应LMS算法的磁悬浮转子谐波振动抑制方法，包括以下步骤：

步骤(1)、建立含质量不平衡和传感器谐波的磁悬浮转子动力学模型

四自由度主被动磁悬浮CMG转子的径向磁轴承控制平动和扭转四自由度，轴向自由度由安装在转子和定子上的永磁环实现无源稳定悬浮。四对径向磁铁和位移传感器对称分布于转子两端A和B，对应4个通道A_x、B_x、A_y、B_y。设转子中心面I，对应设转子A端磁轴承中心面为I₁，B端磁轴承中心面为I₂，转子惯性主轴与I₁,I₂和I分别交于C₁、C₂和C，转子几何轴与I₁,I₂和I分别交于O₁、O₂和O。设两磁轴承中心连线交转子中心面I于磁轴承中心N。

以N为原点在平面I内建立相对于惯性空间的固定坐标系(NXYZ)。在转子中心面I内建立以转子几何中心O点为原点的旋转坐标系Oεη。令l_oc为转子几何中心到转子质心的位移矢量，其中θ为OC与Oε坐标轴的夹角，ψ为OC的长度；相应的，令和分别为由O₁到C₁和由O₂到C₂的矢量，其中φ，为和投影到平面I后与Oε坐标轴的夹角，ξ，ζ为和的模。

主被动磁轴承所提供的轴承力包括两个部分，分为主动磁轴承电磁力和被动磁轴承磁力。A_x通道轴承力f_ax可写为：

f_ax＝f_aex+f_apx

其中，f_aex为A_x通道的主动磁轴承电磁力，f_apx为A_x通道的被动磁轴承磁力。被动磁轴承的磁力大小与位移线性相关，表示为：

f_apx＝K_prx_a

其中，K_pr是被动磁轴承位移刚度，x_a是A_x通道的位移。

当转子在磁中心附近的一定范围内悬浮时，可以将主动磁轴承电磁力线性化为：

f_aex≈K_erx_a+K_ii_ax

其中，K_er、K_i分别为主动磁轴承位移负刚度、电流刚度，i_ax为功放输出电流。

当转子系统含有不平衡质量时，有：

X_a(t)＝x_a(t)+Θ_ax(t)

其中，X_a(t)为转子质心位移，x_a(t)为转子几何中心位移，Θ_ax(t)为质量不平衡引起的位移扰动，记为：

${\mathrm{\Θ}}_{ax}\left(t\right)=l_{o_1{c}_1}cos(\mathrm{\Ω}t+\mathrm{\θ})$

其中，为质量不平衡引起位移扰动的幅值，θ为相位，Ω为转子转速。

在实际系统中，受限于机械加工精度和材料的不均匀特性，磁悬浮转子的位移传感器检测面会出现圆度不理想、材质不均匀、剩磁特性不同等因素，位移传感器的输出信号将会出现多种谐波分量，可表示为：

x_as(t)＝x_a(t)+x_ad(t)

其中，x_ad(t)为传感器谐波，可写为：

$x_{ad}\left(t\right)=\underset{i=1}{\overset{n}{\Σ}}{c}_{ai}sin(i\mathrm{\Ω}t+{\mathrm{\θ}}_i)$

其中，c_ai是传感器谐波分量的幅值，θ_i是传感器谐波分量的相位，n为传感器谐波的最高次数。

将i_ax、X_a、Θ_ax、x_ad依次进行拉普拉斯变换得i_ax(s)、X_a(s)、Θ_ax(s)、x_ad(s)，得到转子动力学方程为：

ms²X_a(s)＝(K_er+K_pr)(X_a(s)-Θ_ax(s))+K_ii_ax(s)

其中，

i_ax(s)＝-K_sG_c(s)G_w(s)(X_a(s)-Θ_ax(s)+x_ad(s))

其中，K_s为位移传感器环节、G_c(s)为控制器环节和G_w(s)为功放环节。则有：

f_ax(s)＝ms²X_a(s)＝(K_er+K_pr)(X_a(s)-Θ_ax(s))-K_sK_iG_c(s)G_w(s)(X_a(s)-Θ_ax(s)+x_ad(s))

通过上式可以得出，等式右边的两项中都包含质量不平衡成分Θ_ax(s)，所以不平衡质量会同时通过控制器和磁轴承产生磁力。而传感器谐波项x_ad(s)仅通过控制器产生电磁力，电磁力中同时包含同频振动力和倍频振动力，而永磁力中只包含同频振动力，所以在进行谐波振动抑制时需要加以区分。

步骤(2)、同频信号中传感器谐波成分的辨识和补偿

由于电磁力和永磁力中都包含同频成分，所以在进行谐波振动抑制时需要分辨同频振动的来源，针对不同的信号源，分别使用不同的方法进行抑制。对于传感器谐波中的同频成分，可以通过控制磁悬浮转子在保护轴承上慢速旋转来近似获取。因为在低转速的状态下，转子质量不平衡所引起的同频位移信号中幅值较小，所以可以认为传感器输出的同频信号均来源于传感器谐波，完成对谐波中同频信号的辨识。在转子高速旋转时，按照慢速旋转下辨识出的同频幅值加入与当前转速相应的同频信号，完成传感器谐波成分的补偿。

步骤(3)、基于频域自适应LMS算法的倍频振动抑制

补偿传感器谐波中的同频信号后，剩下的谐波成分中只含有倍频信号，倍频信号进入控制器，只产生倍频电磁力，因此，抑制倍频振动就要将倍频电流作为控制对象，将传感器谐波引起的倍频振动通过频域自适应LMS算法进行消除，谐波振动抑制的计算过程如下：

使用频域自适应LMS算法进行振动抑制，传感器谐波导致倍频振动的产生，将所需不同频率且与倍频振动信号相关的正弦信号相加作为参考输入，将系统误差作为基本输入信号，在计算过程中，根据相邻两块的融合误差变化情况，改变算法步长和滤波器块长，更好地平衡了收敛速度和稳态精度，在快速收敛的同时，确保得到较小的稳态误差；

步骤(4)、基于陷波器的不平衡振动抑制

补偿传感器谐波中的同频信号，使用频域LMS算法滤除倍频电流，剩下的谐波项为转子的不平衡质量导致的同频电磁力以及同频永磁力，统称为不平衡振动力，对于单一频率信号，可以通过陷波器进行抑制，本发明的控制目标为不平衡振动力，构造不平衡振动输入陷波器，提取同频信号，反馈至控制器从而实现消除不平衡振动力。

本发明基本原理：在控制力矩陀螺的磁悬浮转子系统中，转子质量不平衡和传感器位移误差引起的谐波振动都会大大影响卫星平台的稳定度，灵敏度和分辨率。因此，必须对此进行抑制。首先建立含质量不平衡和传感器谐波的磁悬浮转子动力学模型，通过模型分析谐波振动的主要来源。然后，辨识传感器谐波引起的同频信号，并进行补偿。接着，采用频域自适应LMS算法对传感器谐波引起的倍频振动进行抑制。最后，使用陷波器对质量不平衡引起的同频振动进行抑制。

与目前已有技术相比，本发明的优点在于：提出一种基于频域自适应LMS算法的磁悬浮转子谐波振动抑制方法，通过设置遗忘因子以及计算相邻两块的融合误差，实时进行块长的变化，提高了收敛速度，同时，通过步长的实时调节，减小了系统的稳态误差，有效抑制了磁悬浮转子系统中主要频率成分的谐波振动，适用于存在质量不平衡和传感器谐波的磁悬浮转子谐波振动抑制。

附图说明

图1为本发明的流程图；

图2为A_x通道谐波振动抑制原理框图；

图3为A_y通道谐波振动抑制原理框图；

图4为B_x通道谐波振动抑制原理框图；

图5为B_y通道谐波振动抑制原理框图；

图6为主被动磁悬浮转子系统结构示意图，其中，1为A端磁轴承，2为B端磁轴承，3为转子；

图7为传感器谐波示意图，其中，3为转子，4为传感器；

图8为A_x通道磁轴承控制系统框图；

图9为A_y通道磁轴承控制系统框图；

图10为B_x通道磁轴承控制系统框图；

图11为B_y通道磁轴承控制系统框图；

图12为A_x通道基于频域自适应LMS算法的倍频振动抑制模块；

图13为A_x通道基于陷波器的不平衡振动抑制模块。

具体实施方式

下面结合附图以及具体实施方式进一步说明本发明。

如图1所示，一种基于频域自适应LMS算法的磁悬浮转子谐波振动抑制方法的实施过程是：首先对磁悬浮转子的主要谐波成分进行分析，建立含质量不平衡和传感器谐波的磁悬浮转子动力学模型，然后使用一种频域自适应滤波的方法抑制谐波振动。图2为本发明所采用的磁悬浮转子系统A_x通道谐波振动抑制方法原理框图。图3为本发明所采用的磁悬浮转子系统A_y通道谐波振动抑制方法原理框图。图4为本发明所采用的磁悬浮转子系统B_x通道谐波振动抑制方法原理框图。图5为本发明所采用的磁悬浮转子系统B_y通道谐波振动抑制方法原理框图。谐波振动抑制主要包括传感器同频信号补偿1、倍频振动抑制模块2和不平衡振动抑制模块3。对传感器输出的同频信号进行辨识，将辨识出的传感器同频信号补偿，在此基础上，采用频域LMS算法对倍频振动进行抑制，最后，构造出轴承力，直接以轴承力为控制对象用陷波器抑制不平衡振动。该方法的具体步骤如下：

步骤(1)、建立含质量不平衡和传感器谐波的磁悬浮转子动力学模型

四自由度主被动磁悬浮CMG转子的径向磁轴承控制平动和扭转四自由度，轴向自由度由安装在转子和定子上的永磁环实现无源稳定悬浮。4对径向磁铁和位移传感器对称分布于转子两端A和B，对应4个通道A_x、B_x、A_y、B_y。设转子中心面I，对应设转子A端磁轴承中心面为I₁，B端磁轴承中心面为I₂，转子惯性主轴与I₁,I₂和I分别交于C₁、C₂和C，转子几何轴与I₁,I₂和I分别交于O₁、O₂和O。设两磁轴承中心连线交转子中心面I于磁轴承中心N。

以N为原点在平面I内建立相对于惯性空间的固定坐标系(NXYZ)。在转子中心面I内建立以转子几何中心O点为原点的旋转坐标系Oεη。令l_oc为转子几何中心到转子质心的位移矢量，其中θ为OC与Oε坐标轴的夹角，ψ为OC的长度；相应的，令和分别为由O₁到C₁和由O₂到C₂的矢量，其中φ，为和投影到平面I后与Oε坐标轴的夹角，ξ，ζ为和的模。

主被动磁轴承所提供的轴承力包括两个部分，分为主动磁轴承电磁力和被动磁轴承磁力。以A_x通道为例，轴承力f_ax可写为：

f_ax＝f_aex+f_apx

其中，f_aex为A_x通道的主动磁轴承电磁力，f_apx为A_x通道的被动磁轴承磁力。被动磁轴承的磁力大小与位移线性相关，表示为：

f_apx＝K_prx_a

其中，K_pr是被动磁轴承位移刚度，x_a是A_x通道的位移。

当转子在磁中心附近的一定范围内悬浮时，可以将主动磁轴承电磁力线性化为：

f_aex≈K_erx_a+K_ii_ax

其中，K_er、K_i分别为主动磁轴承位移负刚度、电流刚度，i_ax为功放输出电流。

当转子系统含有不平衡质量时，有：

X_a(t)＝x_a(t)+Θ_ax(t)

其中，X_a(t)为转子质心位移，x_a(t)为转子几何中心位移，Θ_ax(t)为质量不平衡引起的位移扰动，记为：

${\mathrm{\Θ}}_{ax}\left(t\right)=l_{o_1{c}_1}cos(\mathrm{\Ω}t+\mathrm{\θ})$

其中，为质量不平衡引起位移扰动的幅值，θ为相位，Ω为转子转速。

在实际系统中，受限于机械加工精度和材料的不均匀特性，磁悬浮转子的位移传感器检测面会出现圆度不理想、材质不均匀、剩磁特性不同等因素，位移传感器的输出信号将会出现多种谐波分量，可表示为：

x_as(t)＝x_a(t)+x_ad(t)

其中，x_ad(t)为传感器谐波，可写为：

$x_{ad}\left(t\right)=\underset{i=1}{\overset{n}{\Σ}}{c}_{ai}sin(i\mathrm{\Ω}t+{\mathrm{\θ}}_i)$

其中，c_ai是传感器谐波分量的幅值，θ_i是传感器谐波分量的相位，n为传感器谐波的最高次数。

将i_ax、X_a、Θ_ax、x_ad依次进行拉普拉斯变换得i_ax(s)、X_a(s)、Θ_ax(s)、x_ad(s)，得到转子动力学方程为：

ms²X_a(s)＝(K_er+K_pr)(X_a(s)-Θ_ax(s))+K_ii_ax(s)

其中，

i_ax(s)＝-K_sG_c(s)G_w(s)(X_a(s)-Θ_ax(s)+x_ad(s))

其中，K_s为位移传感器环节、G_c(s)为控制器环节和G_w(s)为功放环节。则有：

f_ax(s)＝ms²X_a(s)＝(K_er+K_pr)(X_a(s)-Θ_ax(s))-K_sK_iG_c(s)G_w(s)(X_a(s)-Θ_ax(s)+x_ad(s))

相应的，A_y通道磁轴承控制系统框图如图9所示，A_y通道的轴承力可写为：

f_ay(s)＝ms²Y_a(s)＝(K_er+K_pr)(Y_a(s)-Θ_ay(s))-K_sK_iG_c(s)G_w(s)(Y_a(s)-Θ_ay(s)+y_ad(s))

其中，

B_x通道磁轴承控制系统框图如图10所示，B_x通道的轴承力可写为：

f_bx(s)＝ms²X_b(s)＝(K_er+K_pr)(X_b(s)-Θ_bx(s))-K_sK_iG_c(s)G_w(s)(X_b(s)-Θ_bx(s)+x_bd(s))

其中，

B_y通道磁轴承控制系统框图如图11所示，B_y通道的轴承力可写为：

f_by(s)＝ms²Y_b(s)＝(K_er+K_pr)(Y_b(s)-Θ_by(s))-K_sK_iG_c(s)G_w(s)(Y_b(s)-Θ_by(s)+y_bd(s))

其中，

通过上式可以得出，等式右边的两项中都包含质量不平衡，所以不平衡质量会同时通过控制器和磁轴承产生磁力，而传感器谐波项仅通过控制器产生电磁力，电磁力中同时包含同频振动力和倍频振动力，而永磁力中只包含同频振动力，所以在进行谐波振动抑制时需要加以区分。

步骤(2)、同频信号中传感器谐波成分的辨识和补偿

质量不平衡引起的扰动与转速有关，而传感器谐波引起的同频干扰与转速无关，因此可以利用转速对传感器谐波同频信号辨识并进行补偿。

对于传感器谐波中的同频成分，可以通过控制磁悬浮转子在保护轴承上慢速旋转来近似获取，因为在低转速的状态下，转子质量不平衡所引起的同频位移信号中幅值较小，所以可以认为传感器输出的同频信号均来源于传感器谐波，完成对谐波中同频信号的辨识。在转子高速旋转时，按照慢速旋转下辨识出的同频幅值加入与当前转速相应的同频信号，完成传感器谐波成分的补偿。传感器的实际输出x_as'(t)、y_as'(t)、x_bs'(t)、y_bs'(t)为：

$\left\{\begin{array}{c}{x_{as}}^{\′}\left(t\right)=x_a\left(t\right)+x_{ad}\left(t\right)-{\hat{x}}_{ad}\left(t\right)\\ {y_{as}}^{\′}\left(t\right)=y_a\left(t\right)+y_{ad}\left(t\right)-{\hat{y}}_{ad}\left(t\right)\\ {x_{bs}}^{\′}\left(t\right)=x_b\left(t\right)+x_{bd}\left(t\right)-{\hat{x}}_{bd}\left(t\right)\\ {y_{bs}}^{\′}\left(t\right)=y_b\left(t\right)+y_{bd}\left(t\right)-{\hat{y}}_{bd}\left(t\right)\end{array}\right.$

其中，x_a(t)、y_a(t)为A端X、Y通道传感器的输出；x_b(t)、y_b(t)为B端X、Y通道传感器的输出。为A端X、Y通道的传感器的同频补偿信号， 为A端传感器谐波同频信号幅值和相位的补偿值； 为B端X、Y通道的传感器的同频补偿信号， 为B端传感器谐波同频信号幅值和相位的补偿值。由此可进一步写出补偿直流分量后传感器输出为：

$\left\{\begin{array}{c}{x_{ax}}^{\′}\left(t\right)-x_a\left(t\right)=\underset{i=2}{\overset{n}{\Σ}}{c}_{ai}\sin \left(i\mathrm{\Ω}t+{\mathrm{\θ}}_i\right)+\[c_{a1}\sin \left(\mathrm{\Ω}t+\mathrm{\θ}\right)-{\hat{c}}_{a1}\sin \left(\mathrm{\Ω}t+\widehat{\mathrm{\θ}}\right)\]=\underset{i=2}{\overset{n}{\Σ}}{c}_{ai}\cos \left(i\mathrm{\Ω}t+{\mathrm{\θ}}_i\right)\\ {y_{as}}^{\′}\left(t\right)-y_a\left(t\right)=\underset{i=2}{\overset{n}{\Σ}}{c}_{ai}\cos \left(i\mathrm{\Ω}t+{\mathrm{\θ}}_i\right)+\[c_{a1}\cos \left(\mathrm{\Ω}t+\mathrm{\θ}\right)-{\hat{c}}_{a1}\cos \left(\mathrm{\Ω}t+\widehat{\mathrm{\θ}}\right)\]=\underset{i=2}{\overset{n}{\Σ}}{c}_{ai}\cos \left(i\mathrm{\Ω}t+{\mathrm{\θ}}_i\right)\\ {x_{bs}}^{\′}\left(t\right)-x_b\left(t\right)=\underset{i=2}{\overset{n}{\Σ}}{c}_{bi}\sin \left(i\mathrm{\Ω}t+{\mathrm{\θ}}_i\right)+\[c_{b1}\sin \left(\mathrm{\Ω}t+\mathrm{\θ}\right)-{\hat{c}}_{b1}\sin \left(\mathrm{\Ω}t+\widehat{\mathrm{\θ}}\right)\]=\underset{i=2}{\overset{n}{\Σ}}{c}_{bi}\sin \left(i\mathrm{\Ω}t+{\mathrm{\θ}}_i\right)\\ {y_{bs}}^{\′}\left(t\right)-y_b\left(t\right)=\underset{i=2}{\overset{n}{\Σ}}{c}_{bi}\cos \left(i\mathrm{\Ω}t+{\mathrm{\θ}}_i\right)+\[c_{b1}\cos \left(\mathrm{\Ω}t+\mathrm{\θ}\right)-{\hat{c}}_{b1}\cos \left(\mathrm{\Ω}t+\widehat{\mathrm{\θ}}\right)\]=\underset{i=2}{\overset{n}{\Σ}}{c}_{bi}\cos \left(i\mathrm{\Ω}t+{\mathrm{\θ}}_i\right)\end{array}\right.$

从上式看出，对传感器进行同频信号补偿后，除去直流位移信号，传感器的输出只有倍频信号。

步骤(3)、基于频域自适应LMS算法的倍频振动抑制

补偿传感器谐波中的同频信号后，剩下的谐波成分中只含有倍频信号，倍频信号进入控制器，只产生倍频电磁力，因此，抑制倍频振动就要将倍频电流作为控制对象，将传感器谐波引起的倍频振动通过频域自适应LMS算法进行消除。以A_x通道为例，倍频振动只来源于传感器谐波引起的电磁力，因此本发明以电流为控制目标抑制倍频振动，以电流为基本输入，将所需不同频率且与倍频振动信号相关的正弦信号相加作为参考输入，模块的输出反馈至控制系统的功放输入端，其倍频振动抑制模块如图12所示。

本发明中，所有时域变量均用小写字母表示，所有频域变量均用大写字母表示，向量和矩阵使用黑体表示，标量用斜体表示，傅里叶变换用F(·)表示，傅里叶逆变换用F^-1(·)表示。

设滤波器抽头和块长为N_ax。输入信号通过串并转换变为大小为N_ax的块。定义第k块参考输入u_ax(k)、基本输入e_ax(k)、滤波器输出out_ax(k)、N_ax阶滤波器抽头权向量w_ax(k)依次表示为：

u_ax(k)＝[u_ax(kN_ax-N_ax),u_ax(kN_ax-N_ax+1),…,u_ax(kN_ax+N_ax-1)]^T

e_ax(k)＝[e_ax(kN_ax),e_ax(kN_ax+1),…,e_ax(kN_ax+N_ax-1)]^T

$\begin{array}{c}{\mathrm{out}}_{ax}\left(k\right)={\[{\mathrm{out}}_{ax}\left({\mathrm{kN}}_{ax}-N_{ax}\right),{\mathrm{out}}_{ax}\left({\mathrm{kN}}_a-N_{ax}+1\right),\mathrm{..},{\mathrm{out}}_{ax}\left({\mathrm{kN}}_{ax}+N_{ax}-1\right)\]}^T\\ w_{ax}\left(k\right)={\[w_{0,ax}\left(k\right),w_{1,ax}\left(k\right),\mathrm{...},w_{N_{ax}-1,ax}\left(k\right)\]}^T\end{array}$

滤波器输出out_ax(k)为：

out_ax(k)＝k_axF^-1OUT_ax(k)＝k_axF^-1[U_ax(k)W_ax(k)]

其中，U_ax(k)＝diag{F[u_ax(k)]}，k_ax是N_ax×2N_ax阶约束矩阵：

$k_{ax}=\[0_{N_{ax}},I_{N_{ax}}\]$

其中，是N_ax×N_ax阶零阵，是N_ax×N_ax阶单位阵。

根据随机梯度下降原则，滤波器抽头向量权值更新写为：

W_ax(k+1)＝W_ax(k)+μ_ax(k)F{g_axF^-1[U_ax ^H(k)E_ax(k)]}

其中，W_ax(k+1)为计算出的第k+1块的频域权值向量，U_ax ^H(k)为参考输入频域向量的共轭，E_ax(k)为第k块的频域误差向量，μ_ax(k)为步长，g_ax为梯度约束矩阵：

$g_{ax}=\left[\begin{array}{cc}{I}_{N_{ax}}& 0_{N_{ax}}\\ 0_{N_{ax}}& 0_{N_{ax}}\end{array}\right]$

步长μ_ax(k)控制滤波器抽头权向量从算法的当前迭代到下一次迭代的增量变化，为了保持系统稳定性，需满足：

$0<{\mathrm{\&mu;}}_{ax}\left(k\right)<\frac{1}{N_{ax}{\mathrm{\&lambda;}}_{max}}$

其中，λ_max是参考输入信号相关矩阵R_ax＝E[u_ax(n)u_ax(n)^T]的最大特征值。步长μ_ax(k)的不同取值会影响收敛速度和精度，常采用平均时间常数τ_mse,av和失调系数M这两个量作为衡量频域LMS算法性能的指标。平均时间常数τ_mse,av写为：

${\mathrm{\&tau;}}_{mse,av}=\frac{N_{ax}}{4{\mathrm{\&mu;}}_{ax}\left(k\right){\mathrm{\&lambda;}}_{av}}$

式中，λ_av是矩阵R_ax的特征值的平均值，平均时间常数能反映系统的收敛速度。失调系数M写为：

$M=\frac{{\mathrm{\&mu;}}_{ax}\left(k\right)}{N_{ax}}tr\&lsqb;R_{ax}\&rsqb;$

其中，tr[R_ax]是矩阵R_ax的迹。失调系数可以反映稳态误差。

结合平均时间常数τ_mse,av和失调系数M的表达式可以看出，块长N_ax和步长μ_ax(k)均会影响收敛速度和稳态误差。因此，为了平衡收敛速度和稳态误差这一矛盾，可以实时地改变块长和步长，以达到更快的收敛速度和更小的稳态误差。在初始时刻，可以采用小块长、大步长提高收敛速度，随着误差不断减小，增大块长减小步长来减小稳态误差。具体方法如下：

假设当前块长为N_ax，下一块块长为N_ax'，误差e_ax,ave(k)记为：

$e_{ax,ave}\left(k\right)={\mathrm{\&alpha;e}}_{ax,ave}\left(k-1\right)+\left(1-\mathrm{\&alpha;}\right)\sqrt{\frac{\underset{i=0}{\overset{N_{ax}-1}{\&Sigma;}}{\left(e_{ax}\left({\mathrm{kN}}_{ax}+i\right)\right)}^2}{N_{ax}}}$

其中，e_ax,ave(k)是第k块的融合误差，e_ax,ave(k-1)是第k-1块的融合误差，为第k块的均方误差，第k块的融合误差由第k-1块的融合误差和第k块的均方误差加权得到，α是一个常数且0＜α＜1，e_ax(j)表示第j时刻A_x通道的基本输入，j＝kN_ax,kN_ax+1,…,kN_ax+N_ax-1。块长更新的主要思想是：将当前块的融合误差与上一块的融合误差做比较，如果当前块的融合误差与上一块相比较小，则当前块长增大；如果当前块的融合误差大于上一块，则当前块长减小，但是，由于实际计算中相邻两块的融合误差不可能完全相等，所以每次进行误差比较后块长都会按照更新规则改变，使得计算量大幅增加，因此，在进行前后两块误差的对比时需要留有一定变化裕度，所以，块长的改进更新算法如下：

若β₁e_ax,ave(k-1)＜e_ax,ave(k)＜β₂e_ax,ave(k-1)，则认为误差的变化在允许的范围内，N_ax'＝N_ax。其中，β₁、β₂为常数，且0＜β₁＜1，β₂＞1，两个常数用以保证前后两块误差相比时有一定的裕度，即当前块与上一块的融合误差相差在一定范围内时不用改变块长。

若e_ax,ave(k)≥β₂e_ax,ave(k-1)，则认为当前块的融合误差大于上一块的融合误差，此时，N_ax'＝N_ax/2。

若e_ax,ave(k)≤β₁e_ax,ave(k-1)，则认为当前块的融合误差小于上一块的融合误差，此时，N_ax'＝2N_ax。

块长变化时，权值向量也需要相应地改变。权值向量的更新算法如下：

若N_ax'＝N_ax/2，需要减少N_ax/2个权值向量，有：

${W_{ax}}^{\&prime;}(k+1)={\&lsqb;W_{ax,0}(k+1),W_{ax,2}(k+1),...,W_{ax,2N_{ax}}(k+1)\&rsqb;}^T$

即每隔一个舍弃一个权值向量，W_ax,i(k+1)为计算出的第k+1块频域权值向量中序号为i的值，i＝0,2,4,…,2N_ax。

若N_ax'＝2N_ax，需要增加N_ax个权值向量。由于在时域中补零等效于在频域中插值，因此可以通过对原系数时域补零再频域变换得到新的权系数，即：

w_ax'(k+1)＝[w_ax(k+1),0,…,0]^T

W_ax'(k+1)＝F[w_ax'(k+1)]

其中，w_ax'(k+1)为根据块长更新情况得到的第k+1时刻的时域权值向量，w_ax(k+1)为计算出的第k+1时刻的时域权值向量，W_ax'(k+1)为根据块长更新情况得到的第k+1时刻的频域权值向量。

若N_x'＝N_x，权值向量无需改变。

对于数据块中的每一个信号点采用不同步长以更好地改善收敛性能，通过对每个可调权值赋予不同的步长，频域LMS的收敛速度可以得到改善，步长更新算法为：

${\mathrm{\&mu;}}_{ax,i}\left(k\right)=\frac{{\mathrm{\&mu;}}_0\mathrm{\&beta;}\left(k\right)}{\left|\right|U_{ax,i}\left(k\right)|{|}^2+\mathrm{\&gamma;}}$

其中，μ_ax,i(k)为第k块中第i个步长，μ₀起到控制步长大小的作用，β(k)是控制失调的收敛因子，β(k)＝{1-exp[-υ(k)]}，υ(k)的作用是根据当前块的融合误差e_ax,ave(k)和前一块的融合误差e_ax,ave(k-1)以及常数n和m控制步长变化的形状和速度,U_ax,i(k)为第k块中第i个参考输入的频域信号，||·||表示模值，γ参数是一个较小量，保证步长不会因输入信号能量过低而过大。

同理可得A_y通道的更新算法。滤波器抽头向量权值更新表示为：

W_ay(k+1)＝W_ay(k)+μ_ay(k)F{g_ayF^-1[U_ay ^H(k)E_ay(k)]}

其中，W_ay(k)为A_y通道第k块频域权值向量，W_ay(k+1)为计算出的A_y通道第k+1块频域权值向量，U_ay ^H(k)为参考输入频域向量的共轭，U_ay(k)＝diag{F[u_ay(k)]}，u_ay(k)为参考输入的时域向量，E_ay(k)为第k块频域误差向量，μ_ay(k)为步长，g_ay为梯度约束矩阵。滤波器输出out_ay(k)为：

out_ay(k)＝k_ayF^-1OUT_ay(k)＝k_ayF^-1[U_ay(k)W_ay(k)]

其中，OUT_ay(k)为第k块滤波器的频域输出。k_ay是N_ay×2N_ay阶约束矩阵：

$k_{ay}=\&lsqb;0_{N_{ay}},I_{N_{ay}}\&rsqb;$

假设当前块长为N_ay，下一块块长为N_ay'，误差e_ay,ave(k)记为：

$e_{ay,ave}\left(k\right)={\mathrm{\&alpha;e}}_{ay,ave}\left(k-1\right)+\left(1-\mathrm{\&alpha;}\right)\sqrt{\frac{\underset{i=0}{\overset{N_{ay}-1}{\&Sigma;}}{\left(e_{ay}\left({\mathrm{kN}}_{ay}+i\right)\right)}^2}{N_{ay}}}$

其中，e_ay,ave(k)是第k块的融合误差，e_ay,ave(k-1)是第k-1块的融合误差，为第k块的均方误差，第k块的融合误差由第k-1块的融合误差和第k块的均方误差加权得到，α是一个常数且0＜α＜1，e_ay(j)表示第j时刻A_y通道的基本输入，j＝kN_ay,kN_ay+1,…,kN_ay+N_ay-1。块长的具体更新算法如下：

若β₁e_ay,ave(k-1)＜e_ay,ave(k)＜β₂e_ay,ave(k-1)，则认为误差的变化在允许的范围内，N_ay'＝N_ay。

若e_ay,ave(k)≥β₂e_ay,ave(k-1)，则视为当前块的融合误差大于上一块的融合误差，此时，N_ay'＝N_ay/2。

若e_ay,ave(k)≤β₁e_ay,ave(k-1)，则视为当前块的融合误差小于上一块的融合误差，此时，N_ay'＝2N_ay。

权值向量的更新算法如下：

若N_ay'＝N_ay/2，需要减少N_ay/2个权值向量，有：

${W_{ay}}^{\&prime;}(k+1)={\&lsqb;W_{ay,0}(k+1),W_{ay,2}(k+1),...,W_{ay,2N_{ay}}(k+1)\&rsqb;}^T$

即每隔一个舍弃一个权值向量，W_ay,i(k+1)为计算出的第k+1块频域权值向量中序号为i的值，i＝0,2,4,…,2N_ay。

若N_ay'＝2N_ay，需要增加N_ay个权值向量。由于在时域中补零等效于在频域中插值，因此可以通过对原系数时域补零再频域变换得到新的权系数，即：

w_ay'(k+1)＝[w_ay(k+1),0,…,0]^T

W_ay'(k+1)＝F[w_ay'(k+1)]

其中，w_ay'(k+1)为根据块长更新情况得到的第k+1时刻的时域权值向量，w_ay(k+1)为计算出的第k+1时刻的时域权值向量，W_ay'(k+1)为根据块长更新情况得到的第k+1时刻的频域权值向量。

若N_ay'＝N_ay，权值向量无需改变。

步长更新算法为：

${\mathrm{\&mu;}}_{ay,i}\left(k\right)=\frac{{\mathrm{\&mu;}}_0\mathrm{\&beta;}\left(k\right)}{\left|\right|U_{ay,i}\left(k\right)|{|}^2+\mathrm{\&gamma;}}$

其中,

β(k)＝{1-exp[-υ(k)]}

$\mathrm{\&upsi;}\left(k\right)=\left|\right|\frac{e_{ay,ave}\left(k\right)}{e_{ay,ave}(k-1)}|{|}^n.||e_{ay,ave}(k-1)-e_{ay,ave}\left(k\right)|{|}^m$

U_ay,i(k)为第k块中第i个参考输入的频域信号。

同理可得B_x通道的更新算法。滤波器抽头向量权值更新表示为：

W_bx(k+1)＝W_bx(k)+μ_bx(k)F{g_bxF^-1[U_bx ^H(k)E_bx(k)]}

其中，W_bx(k)为B_x通道第k块频域权值向量，W_bx(k+1)为计算出的B_x通道第k+1块频域权值向量，U_bx ^H(k)为参考输入频域向量的共轭，U_bx(k)＝diag{F[u_bx(k)]}，u_bx(k)为参考输入的时域向量，E_bx(k)为第k块频域误差向量，μ_bx(k)为步长，g_bx为梯度约束矩阵。滤波器输出out_bx(k)为：

out_bx(k)＝k_bxF^-1OUT_bx(k)＝k_bxF^-1[U_bx(k)W_bx(k)]

其中，OUT_bx(k)为第k块滤波器的频域输出。k_bx是N_bx×2N_bx阶约束矩阵：

$k_{bx}=\&lsqb;0_{N_{bx}},I_{N_{bx}}\&rsqb;$

假设当前块长为N_bx，下一块块长为N_bx'，误差e_bx,ave(k)记为：

$e_{bx,ave}\left(k\right)={\mathrm{\&alpha;e}}_{bx,ave}(k-1)+(1-\mathrm{\&alpha;})\sqrt{\frac{\underset{i=0}{\overset{N_{bx}-1}{\&Sigma;}}{\left(e_{bx}\right({\mathrm{kN}}_{bx}+i\left)\right)}^2}{N_{bx}}}$

其中，e_bx,ave(k)是第k块的融合误差，e_bx,ave(k-1)是第k-1块的融合误差，为第k块的均方误差，第k块的融合误差由第k-1块的融合误差和第k块的均方误差加权得到，α是一个常数且0＜α＜1，e_bx(j)表示第j时刻B_x通道的基本输入，j＝kN_bx,kN_bx+1,…,kN_bx+N_bx-1。块长的具体更新算法如下：

若β₁e_bx,ave(k-1)＜e_bx,ave(k)＜β₂e_bx,ave(k-1)，则认为误差的变化在允许的范围内，N_bx'＝N_bx。

若e_bx,ave(k)≥β₂e_bx,ave(k-1)，则视为当前块的融合误差大于上一块的融合误差，此时，N_bx'＝N_bx/2。

若e_bx,ave(k)≤β₁e_bx,ave(k-1)，则视为当前块的融合误差小于上一块的融合误差，此时，N_bx'＝2N_bx。

权值向量的更新算法如下：

若N_bx'＝N_bx/2，需要减少N_bx/2个权值向量，有：

${W_{bx}}^{\&prime;}(k+1)={\&lsqb;W_{bx,0}(k+1),W_{bx,2}(k+1),...,W_{bx,2N_{bx}}(k+1)\&rsqb;}^T$

即每隔一个舍弃一个权值向量，W_bx,i(k+1)为计算出的第k+1块频域权值向量中序号为i的值，i＝0,2,4,…,2N_bx。

若N_bx'＝2N_bx，需要增加N_bx个权值向量。由于在时域中补零等效于在频域中插值，因此可以通过对原系数时域补零再频域变换得到新的权系数，即：

w_bx'(k+1)＝[w_bx(k+1),0,…,0]^T

W_bx'(k+1)＝F[w_bx'(k+1)]

其中，w_bx'(k+1)为根据块长更新情况得到的第k+1时刻的时域权值向量，w_bx(k+1)为计算出的第k+1时刻的时域权值向量，W_bx'(k+1)为根据块长更新情况得到的第k+1时刻的频域权值向量。

若N_bx'＝N_bx，权值向量无需改变。

步长更新算法为：

${\mathrm{\&mu;}}_{bx,i}\left(k\right)=\frac{{\mathrm{\&mu;}}_0\mathrm{\&beta;}\left(k\right)}{\left|\right|U_{bx,i}\left(k\right)|{|}^2+\mathrm{\&gamma;}}$

其中,

β(k)＝{1-exp[-υ(k)]}

$\mathrm{\&upsi;}\left(k\right)=\left|\right|\frac{e_{bx,ave}\left(k\right)}{e_{bx,ave}(k-1)}|{|}^n\&CenterDot;||e_{bx,ave}(k-1)-e_{bx,ave}\left(k\right)|{|}^m$

U_bx,i(k)为第k块中第i个参考输入的频域信号。

同理可得B_y通道的更新算法。滤波器抽头向量权值更新表示为：

W_by(k+1)＝W_by(k)+μ_by(k)F{g_byF^-1[U_by ^H(k)E_by(k)]}

其中，W_by(k)为B_y通道第k块频域权值向量，W_by(k+1)为计算出的B_y通道第k+1块频域权值向量，U_by ^H(k)为参考输入频域向量的共轭，U_by(k)＝diag{F[u_by(k)]}，u_by(k)为参考输入的时域向量，E_by(k)为第k块频域误差向量，μ_by(k)为步长，g_by为梯度约束矩阵。滤波器输出out_by(k)为：

out_by(k)＝k_byF^-1OUT_by(k)＝k_byF^-1[U_by(k)W_by(k)]

其中，OUT_by(k)为第k块滤波器的频域输出。k_by是N_by×2N_by阶约束矩阵：

$k_{by}=\&lsqb;0_{N_{by}},I_{N_{by}}\&rsqb;$

假设当前块长为N_by，下一块块长为N_by'，误差e_by,ave(k)记为：

$e_{by,ave}\left(k\right)={\mathrm{\&alpha;e}}_{by,ave}(k-1)+(1-\mathrm{\&alpha;})\sqrt{\frac{\underset{i=0}{\overset{N_{by}-1}{\&Sigma;}}{\left(e_{by}\right({\mathrm{kN}}_{by}+i\left)\right)}^2}{N_{by}}}$

其中，e_by,ave(k)是第k块的融合误差，e_by,ave(k-1)是第k-1块的融合误差，为第k块的均方误差，第k块的融合误差由第k-1块的融合误差和第k块的均方误差加权得到，α是一个常数且0＜α＜1，e_by(j)表示第j时刻B_y通道的基本输入，j＝kN_by,kN_by+1,…,kN_by+N_by-1。块长的具体更新算法如下：

若β₁e_by,ave(k-1)＜e_by,ave(k)＜β₂e_by,ave(k-1)，则认为误差的变化在允许的范围内，N_by'＝N_by。

若e_by,ave(k)≥β₂e_by,ave(k-1)，则视为当前块的融合误差大于上一块的融合误差，此时，N_by'＝N_by/2。

若e_by,ave(k)≤β₁e_by,ave(k-1)，则视为当前块的融合误差小于上一块的融合误差，此时，N_by'＝2N_by。

权值向量的更新算法如下：

若N_by'＝N_by/2，需要减少N_by/2个权值向量，有：

${W_{by}}^{\&prime;}(k+1)={\&lsqb;W_{by,0}(k+1),W_{by,2}(k+1),...,W_{by,2N_{by}}(k+1)\&rsqb;}^T$

即每隔一个舍弃一个权值向量，W_by,i(k+1)为计算出的第k+1块频域权值向量中序号为i的值，i＝0,2,4,…,2N_by。

若N_by'＝2N_by，需要增加N_by个权值向量。由于在时域中补零等效于在频域中插值，因此可以通过对原系数时域补零再频域变换得到新的权系数，即：

w_by'(k+1)＝[w_by(k+1),0,…,0]^T

W_by'(k+1)＝F[w_by'(k+1)]

其中，w_by'(k+1)为根据块长更新情况得到的第k+1时刻的时域权值向量，w_by(k+1)为计算出的第k+1时刻的时域权值向量，W_by'(k+1)为根据块长更新情况得到的第k+1时刻的频域权值向量。

若N_by'＝N_by，权值向量无需改变。

步长更新算法为：

${\mathrm{\&mu;}}_{by,i}\left(k\right)=\frac{{\mathrm{\&mu;}}_0\mathrm{\&beta;}\left(k\right)}{\left|\right|U_{by,i}\left(k\right)|{|}^2+\mathrm{\&gamma;}}$

其中,

β(k)＝{1-exp[-υ(k)]}

$\mathrm{\&upsi;}\left(k\right)=\left|\right|\frac{e_{by,ave}\left(k\right)}{e_{by,ave}\left(k-1\right)}|{|}^n\&CenterDot;||e_{by,ave}\left(k-1\right)-e_{by,ave}\left(k\right)|{|}^m$

U_by,i(k)为第k块中第i个参考输入的频域信号。

步骤(4)、基于陷波器的不平衡振动抑制

补偿传感器谐波中的同频信号，使用频域LMS算法滤除倍频电流，剩下的谐波项为转子的不平衡质量导致的同频电磁力以及同频永磁力，统称为不平衡振动力，对于单一频率信号，可以通过陷波器进行抑制，本发明的控制目标为不平衡振动力，以A_x通道为例，将构造出的不平衡振动力F_ax(s)作为陷波器的输入，输出反馈至控制器的输入。构造出的不平衡振动力F_ax(s)可写为：

F_ax(s)＝K_ii_ax(s)+(K_er+K_pr)x_a(s)

如图13所示，跟踪滤波器N_f(s)的传递函数为：

$N_f\left(s\right)=\frac{s}{s^2+{\mathrm{\&Omega;}}^2}$

以质量不平衡Θ_ax(s)为输入，轴承力F_ax(s)为输出，对应的传递函数为：

$\frac{F_{ax}\left(s\right)}{{\mathrm{\&Theta;}}_{ax}\left(s\right)}=\frac{K_{er}+K_{pr}-K_i{K}_s{G}_c\left(s\right)G_w\left(s\right)}{1+K_i{K}_s{G}_c\left(s\right)G_w\left(s\right)P\left(s\right)-(K_{er}+K_{pr})P\left(s\right)+{\mathrm{\&epsiv;K}}_i{N}_f\left(s\right)G_w\left(s\right)}$

其中，P(s)为转子系统传递函数。则，有：

$\underset{s\&RightArrow;j\mathrm{\&Omega;}}{\lim }{F}_{ax}\left(j\mathrm{\&Omega;}\right)=0$

由上式可得，该陷波器可以对不平衡振动进行抑制。

同理可得A_y通道不平衡振动抑制方法。将构造出的不平衡振动力F_ay(s)作为陷波器的输入，输出反馈至控制器的输入。构造出的不平衡振动力F_ay(s)可写为：

F_ay(s)＝K_ii_ay(s)+(K_er+K_pr)y_a(s)

以质量不平衡Θ_ay(s)为输入，轴承力F_ay(s)为输出，对应的传递函数为：

$\frac{F_{ay}\left(s\right)}{{\mathrm{\&Theta;}}_{ay}\left(s\right)}=\frac{K_{er}+K_{pr}-K_i{K}_s{G}_c\left(s\right)G_w\left(s\right)}{1+K_i{K}_s{G}_c\left(s\right)G_w\left(s\right)P\left(s\right)-(K_{er}+K_{pr})P\left(s\right)+{\mathrm{\&epsiv;K}}_i{N}_f\left(s\right)G_w\left(s\right)}$

则，有：

$\underset{s\&RightArrow;j\mathrm{\&Omega;}}{\lim }{F}_{ay}\left(j\mathrm{\&Omega;}\right)=0$

由上式可得，该陷波器可以对不平衡振动进行抑制。

同理可得B_x通道不平衡振动抑制方法。将构造出的不平衡振动力F_bx(s)作为陷波器的输入，输出反馈至控制器的输入。构造出的不平衡振动力F_bx(s)可写为：

F_bx(s)＝K_ii_bx(s)+(K_er+K_pr)x_b(s)

以质量不平衡Θ_bx(s)为输入，轴承力F_y(s)为输出，对应的传递函数为：

$\frac{F_{bx}\left(s\right)}{{\mathrm{\&Theta;}}_{bx}\left(s\right)}=\frac{K_{er}+K_{pr}-K_i{K}_s{G}_c\left(s\right)G_w\left(s\right)}{1+K_i{K}_s{G}_c\left(s\right)G_w\left(s\right)P\left(s\right)-\left(K_{er}+K_{pr}\right)P\left(s\right)+{\mathrm{\&epsiv;K}}_i{N}_f\left(s\right)G_w\left(s\right)}$

则，有：

$\underset{s\&RightArrow;j\mathrm{\&Omega;}}{\lim }{F}_{bx}\left(j\mathrm{\&Omega;}\right)=0$

由上式可得，该陷波器可以对不平衡振动进行抑制。

同理可得B_y通道不平衡振动抑制方法。将构造出的不平衡振动力F_by(s)作为陷波器的输入，输出反馈至控制器的输入。构造出的不平衡振动力F_by(s)可写为：

F_by(s)＝K_ii_by(s)+(K_er+K_pr)y_b(s)

以质量不平衡Θ_by(s)为输入，轴承力F_y(s)为输出，对应的传递函数为：

$\frac{F_{by}\left(s\right)}{{\mathrm{\&Theta;}}_{by}\left(s\right)}=\frac{K_{er}+K_{pr}-K_i{K}_s{G}_c\left(s\right)G_w\left(s\right)}{1+K_i{K}_s{G}_c\left(s\right)G_w\left(s\right)P\left(s\right)-(K_{er}+K_{pr})P\left(s\right)+{\mathrm{\&epsiv;K}}_i{N}_f\left(s\right)G_w\left(s\right)}$

则，有：

$\underset{s\&RightArrow;j\mathrm{\&Omega;}}{\lim }{F}_{by}\left(j\mathrm{\&Omega;}\right)=0$

由上式可得，该陷波器可以对不平衡振动进行抑制。

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

Claims (2)

1.一种基于频域自适应LMS算法的磁悬浮转子谐波振动抑制方法，其特征在于：该方法包括以下步骤：

步骤(1)、建立包含质量不平衡和传感器谐波的磁悬浮转子动力学模型

四自由度主被动磁悬浮CMG(ControlMomentGyro)转子的径向磁轴承控制平动和扭转四个自由度，轴向自由度由安装在转子和定子上的永磁环实现无源稳定悬浮，4对径向磁铁和位移传感器对称分布于转子两端A和B，对应4个通道A_x、B_x、A_y、B_y，设转子中心面I，对应设转子A端磁轴承中心面为I₁，B端磁轴承中心面为I₂，转子惯性主轴与I₁,I₂和I分别交于C₁、C₂和C，转子几何轴与I₁,I₂和I分别交于O₁、O₂和O，设两磁轴承中心连线交转子中心面I于磁轴承中心N；

以N为原点在平面I内建立相对于惯性空间的固定坐标系(NXYZ)，在转子中心面I内建立以转子几何中心O点为原点的旋转坐标系Oεη，令l_oc为转子几何中心到转子质心的位移矢量，其中θ为OC与Oε坐标轴的夹角，ψ为OC的长度；相应的，令和分别为由O₁到C₁和由O₂到C₂的矢量，其中φ，为和投影到平面I后与Oε坐标轴的夹角，ξ，ζ为和的模；

磁轴承所提供的轴承力包括两个部分，分为主动磁轴承电磁力和被动磁轴承磁力，A_x通道的轴承力f_ax可写为：

f_ax＝f_aex+f_apx

其中，f_aex为A_x通道的主动磁轴承电磁力，f_apx为A_x通道的被动磁轴承磁力，被动磁轴承的磁力大小与位移线性相关，表示为：

f_apx＝K_prx_a

其中，K_pr是被动磁轴承位移刚度，x_a是A_x通道的位移；

当转子在磁中心附近的一定范围内悬浮时，可以将主动磁轴承电磁力线性化为：

f_aex≈K_erx_a+K_ii_ax

其中，K_er、K_i分别为主动磁轴承位移负刚度、电流刚度，i_ax为功放输出电流；

当转子系统含有不平衡质量时，有：

X_a(t)＝x_a(t)+Θ_ax(t)

其中，X_a(t)为转子质心位移，x_a(t)为转子几何中心位移，Θ_ax(t)为质量不平衡引起的位移扰动，记为：

${\mathrm{\&Theta;}}_{ax}\left(t\right)=l_{o_1{c}_1}cos(\mathrm{\&Omega;}t+\mathrm{\&theta;})$

其中，为质量不平衡引起位移扰动的幅值，θ为相位，Ω为转子转速；

在实际系统中，受限于机械加工精度和材料的不均匀特性，磁悬浮转子的位移传感器检测面会出现圆度不理想、材质不均匀、剩磁特性不同等因素，位移传感器的输出信号将会出现多种谐波分量，可表示为：

x_as(t)＝x_a(t)+x_ad(t)

其中，x_ad(t)为传感器谐波，可写为：

$x_{ad}\left(t\right)=\underset{i=1}{\overset{n}{\&Sigma;}}{c}_{ai}sin(i\mathrm{\&Omega;}t+{\mathrm{\&theta;}}_i)$

其中，c_ai是传感器谐波分量的幅值，θ_i是传感器谐波分量的相位，n为传感器谐波的最高次数；

将i_ax、X_a、Θ_ax、x_ad依次进行拉普拉斯变换得i_ax(s)、X_a(s)、Θ_ax(s)、x_ad(s)，得到转子动力学方程为：

ms²X_a(s)＝(K_er+K_pr)(X_a(s)-Θ_ax(s))+K_ii_ax(s)

其中，

i_ax(s)＝-K_sG_c(s)G_w(s)(X_a(s)-Θ_ax(s)+x_ad(s))

其中，K_s为位移传感器环节、G_c(s)为控制器环节和G_w(s)为功放环节，则有：

f_ax(s)＝ms²X_a(s)＝(K_er+K_pr)(X_a(s)-Θ_ax(s))-K_sK_iG_c(s)G_w(s)(X_a(s)-Θ_ax(s)+x_ad(s))

通过上式可以得出，等式右边的两项中都包含质量不平衡成分Θ_ax(s），所以不平衡质量会同时通过控制器和磁轴承产生磁力，而传感器谐波项x_ad(s)仅通过控制器产生电磁力，电磁力中同时包含同频振动力和倍频振动力，而永磁力中只包含同频振动力，所以在进行谐波振动抑制时需要加以区分；

步骤(2)、同频信号中传感器谐波成分的辨识和补偿

由于电磁力和永磁力中都包含同频成分，所以在进行谐波振动抑制时需要分辨同频振动的来源，针对不同的信号源，分别使用不同的方法进行抑制，对于传感器谐波中的同频成分，可以通过控制磁悬浮转子在保护轴承上慢速旋转来近似获取，因为在低转速的状态下，转子质量不平衡所引起的同频位移信号中幅值较小，所以可以认为传感器输出的同频信号均来源于传感器谐波，完成对谐波中同频信号的辨识，在转子高速旋转时，按照慢速旋转下辨识出的同频幅值加入与当前转速相应的同频信号，完成传感器谐波成分的补偿；

步骤(3)、基于频域自适应LMS算法的倍频振动抑制

补偿传感器谐波中的同频信号后，剩下的谐波成分中只含有倍频信号，倍频信号进入控制器，只产生倍频电磁力，因此，抑制倍频振动就要将倍频电流作为控制对象，将传感器谐波引起的倍频振动通过频域自适应LMS(LeastMeanSquare)算法进行消除，将所需不同频率且与倍频振动信号相关的正弦信号相加作为参考输入，将系统误差作为基本输入信号，在计算过程中，根据相邻两块的融合误差的变化情况，改变算法步长和滤波器块长，更好地平衡了收敛速度和稳态精度，在快速收敛的同时，确保得到较小的稳态误差；

步骤(4)、基于陷波器的不平衡振动抑制

补偿传感器谐波中的同频信号，使用频域LMS算法滤除倍频电流，剩下的谐波项为转子的不平衡质量导致的同频电磁力以及同频永磁力，统称为不平衡振动力，对于单一频率信号，可以通过陷波器进行抑制，本发明的控制目标为不平衡振动力，构造不平衡振动输入陷波器，提取同频信号，反馈至控制器从而实现消除不平衡振动力。

2.根据权利要求1所述的一种基于频域自适应LMS算法的磁悬浮转子谐波振动抑制方法，其特征在于：所述的步骤(3)块长更新算法为：

假设当前块长为N_ax，下一块块长为N_ax'，融合误差e_ax,ave(k)记为：

$e_{ax,ave}\left(k\right)={\mathrm{\&alpha;e}}_{ax,ave}(k-1)+(1-\mathrm{\&alpha;})\sqrt{\frac{\underset{i=0}{\overset{N_{ax}-1}{\&Sigma;}}{\left(e_{ax}\right({\mathrm{kN}}_{ax}+i\left)\right)}^2}{N_{ax}}}$

其中，e_ax,ave(k)是第k块的融合误差，e_ax,ave(k-1)是第k-1块的融合误差，为第k块的均方误差，第k块的融合误差由第k-1块的融合误差和第k块的均方误差加权得到，α是一个常数且0＜α＜1，e_ax(j)表示第j时刻A_x通道的基本输入，j＝kN_ax,kN_ax+1,…,kN_ax+N_ax-1，块长更新的主要思想是：将当前块的融合误差与上一块的融合误差做比较，如果当前块的融合误差与上一块相比较小，则当前块长增大；如果当前块的融合误差大于上一块，则当前块长减小，但是，由于实际计算中相邻两块的融合误差不可能完全相等，所以每次进行误差比较后块长都会按照更新规则改变，使得计算量大幅增加，因此，在进行前后两块误差的对比时需要留有一定变化裕度，所以，块长的改进更新算法如下：

若β₁e_ax,ave(k-1)＜e_ax,ave(k)＜β₂e_ax,ave(k-1)，则认为误差的变化在允许的范围内，块长不变，N_ax'＝N_ax，其中，β₁、β₂为两个常数，0＜β₁＜1且β₂＞1，两个常数用以保证前后两块误差相比时有一定的裕度，即当前块与上一块的融合误差相差在一定范围内时不用改变块长；

若e_ax,ave(k)≥β₂e_ax,ave(k-1)，则认为当前块的融合误差大于上一块的融合误差，此时，N_ax'＝N_ax/2；

若e_ax,ave(k)≤β₁e_ax,ave(k-1)，则认为当前块的融合误差小于上一块的融合误差，此时，N_ax'＝2N_ax；

所述的步骤(3)步长更新算法为：

对于数据块中的每一个信号点采用不同步长以更好地改善收敛性能，通过对每个可调权值赋予不同的步长，频域LMS的收敛速度可以得到改善，步长更新算法为：

${\mathrm{\&mu;}}_{ax,i}\left(k\right)=\frac{{\mathrm{\&mu;}}_0\mathrm{\&beta;}\left(k\right)}{\left|\right|U_{ax,i}\left(k\right)|{|}^2+\mathrm{\&gamma;}}$

其中，μ_ax,i(k)为第k块中第i个步长，μ₀起到控制步长大小的作用，β(k)是控制失调的收敛因子，β(k)＝{1-exp[-υ(k)]}，υ(k)的作用是根据当前块的融合误差e_ax,ave(k)和前一块的融合误差e_ax,ave(k-1)以及常数n和m控制步长变化的形状和速度，U_ax,i(k)为第k块中第i个参考输入的频域信号，||·||表示模值，γ参数是一个较小量，保证步长不会因输入信号能量过低而过大。