一种基于LQR算法的旋翼振动主动控制方法

摘要

本发明公开了一种基于LQR算法的旋翼振动主动控制方法，该方法首先对旋翼模型进行有限元的模态分析，求出刚度矩阵、质量矩阵和阻尼矩阵，然后将控制方程写成状态空间的表示形式。利用最优控制理论中的LQR算法，把旋翼的振动主动控制问题等价成输出调节器问题。采用基于LQR算法的振动主动控制之后，旋翼的振动剧烈程度明显改善。本发明是主动控制理论与有限元分析的联合应用，便于在旋翼的设计阶段进行振动控制预估，提高了分析效率，通过采用主动控制的方法，改善了旋翼的振动特性，对于具体应用中有很强工程实践意义。

说明书

技术领域

本发明适用于旋翼振动的主动控制领域。旋翼在使用实践中面临振动过于剧烈而导致减少使用寿命、振动导致噪声等问题。由于直升机在通用航空中的使用愈发广泛，旋翼剧烈振动引发的问题愈发突出。本发明通过采用施加主动控制力矩的方式，结合控制理论中的LQR算法，与有限元方法联合，达到控制旋翼振动得目的。

背景技术

本发明公开了一种基于LQR算法的旋翼振动主动控制方法。要求旋翼外表面铺设有压电纤维铺层，用以在通电情况下产生控制力矩。该方法首先对旋翼模型进行有限元的模态分析，求出刚度矩阵、质量矩阵和阻尼矩阵，然后将控制方程写成状态空间的表示形式。利用控制理论中的LQR算法，把旋翼的振动主动控制问题等价成输出调节器问题。采用基于LQR算法的振动控制之后，旋翼的振动剧烈程度明显改善。本发明是主动控制理论与有限元分析的联合应用，便于在旋翼的设计阶段进行振动控制预估，提高了分析效率，通过采用主动施加力矩的控制方法，改善了旋翼的振动特性，对于具体应用中具有很强工程实践意义。

发明内容

本发明的技术解决问题：利用控制理论中的LQR算法，把旋翼的振动主动控制问题等价成输出调节器问题。采用施加控制力矩的主动控制方式，结合控制理论中的LQR算法，与有限元分析联合，达到控制旋翼振动特性的目的。

本发明技术解决方案：首先对旋翼模型进行有限元的模态分析，求出刚度矩阵、质量矩阵和阻尼矩阵，然后将控制方程写成状态空间的表示形式，这样就可以利用现代控制理论中的先进控制方法。利用控制理论中的LQR算法，把旋翼的振动主动控制问题等价成输出调节器问题。通过LQR方法计算得到状态反馈增益向量，判断进行控制之后的旋翼振动特性是否满足约束条件。如果不满足，重新确定控制力矩的个数以及各个控制力矩的施加位置，重复LQR算法，直到旋翼的振动特性满足约束条件为止。

本发明为一种基于LQR算法的旋翼振动主动控制方法，实现步骤包括：

第一步：首先对旋翼模型进行有限元的模态分析，得到旋翼的刚度矩阵K₀，质量矩阵M₀，以及阻尼矩阵C₀；以及有限元模型中的全部位移坐标列向量x；

第二步，在旋翼外表面铺设压电纤维铺层，对压电纤维铺层通电产生控制力矩，通过对旋翼施加控制力矩的方式达到控制振动特性的目的；

第三步，确定控制力矩的个数r以及各个控制力矩的施加位置，所有控制力矩中的力的大小组成列向量u；

第四步：根据第一步得到的刚度矩阵K₀、质量矩阵M₀，阻尼矩阵C₀，以及全部位移坐标列向量x和第三步确定的控制力矩个数和位置，得到旋翼的振动控制微分方程；将旋翼振动控制微分方程改写成标准的状态空间方程表示形式；第三步中施加控制力矩位置处的节点位移作为输出响应向量y；

第五步：根据第四步得到的状态空间方程，采用控制理论中的LQR算法，并且基于输出调节器的最优控制理论，进行旋翼振动的主动控制设计；

第六步：通过第五步中的LQR算法计算得到状态反馈增益矩阵G；通过G可以求出振动过程中施加的控制力矩的大小，从而开展旋翼振动主动控制；

第七步：基于第六步中旋翼实施振动主动控制之后的振动特性效果，判断进行主动控制之后的旋翼振动特性是否满足约束条件；如果满足，则结束；如果不满足，返回第三步，控制力矩的个数由原来的r变为r+1，然后重复第四步到第六步操作，直到旋翼的振动特性满足约束条件为止。

所述第三步中，r个力矩的位置分别位于旋翼展向长度的r+1等分的r个等分点处，对于r＝2时，力矩位置分别位于将旋翼长度进行3等分的两个等分点处。

所述第五步和第六步中，采用控制理论中的LQR算法，并且基于输出调节器的最优控制理论，进行旋翼振动的主动控制，及通过LQR方法计算得到状态反馈增益矩阵G的具体步骤如下：

(1)在输出调节器问题中，性能指标取为二次型形式如下：

$J=\frac{1}{2}{\int }_0^{\infty }[y^{T}Qy+u^{T}Ru]dt$

其中，利用第三步所有控制力矩中的力的大小组成列向量u和第四步输出响应向量y。最优控制理论要求J取最小值时的控制是最优控制。

控制理论中的LQR算法，要求矩阵Q和矩阵R是正定常数的权系数矩阵。计算中均取为单位对角矩阵。

(2)根据矩阵Q和矩阵R以及第四步中的状态空间方程，解矩阵代数的Riccati方程，可以计算得到状态反馈增益矩阵G。(Riccati方程是控制理论的普遍方程)

本发明与现有技术相比的优点在于：有限元方法与控制理论相结合，充分发挥多学科综合的优势，避免了单一有限元分析的盲目性以及单一控制理论中辨识技术的不准确性。在有限元分析应用的基础之上，采用基于LQR算法进行振动控制，旋翼的振动剧烈程度明显改善，从而为工程应用中的振动特性设计提供指导，不但节省分析成本，也提高了控制效率。

附图说明

图1是本发明基于LQR算法的旋翼振动主动控制的实施流程；

图2是本发明中旋翼有限元模型的示意图；

图3是本发明中旋翼某一横截面的示意图，外部蒙皮处布置有压电纤维铺层；

图4是本发明实施例中施加控制力矩的示意图，共计3个控制力矩；

图5是本发明实施例中无振动控制时的旋翼端点位移响应；

图6是本发明实施例中有LQR振动主动控制时的旋翼端点位移响应；

图7是本发明实施例中有无控制时的旋翼端点位移响应的对比。

具体实施方式

如图1所示，本发明具体实现步骤是：

第一步：首先对旋翼模型进行有限元的模态分析，得到旋翼的刚度矩阵K₀，质量矩阵M₀，以及阻尼矩阵C₀；以及有限元模型中的全部位移坐标列向量x；

有限元分析采用一般的商业有限元软件即可。这一步需要对有限元模态分析中的刚度矩阵K₀，质量矩阵M₀，阻尼矩阵C₀还有全部位移坐标列向量x进行输出操作。

第二步：在旋翼外表面铺设压电纤维铺层，对压电纤维铺层通电产生控制力矩，通过对旋翼施加控制力矩的方式达到控制振动特性的目的；

第三步：确定控制力矩的个数r以及各个控制力矩的施加位置，所有控制力矩中的力的大小组成列向量u；其中r个力矩的位置分别位于旋翼展向长度的r+1等分的r个等分点处，对于r＝2时，力矩位置分别位于将旋翼长度进行3等分的两个等分点处。

因为压电纤维铺层在通电状态下产生的力是内力，不是外力。不同位置处的压电纤维铺层在通电状态下对旋翼结构可以形成力矩的作用效果。力矩可以等效成对于结构的大小相等、方向相反、不在同一处作用点的一对作用力的作用效果。

第四步：根据第一步得到的刚度矩阵K₀、质量矩阵M₀，阻尼矩阵C₀，以及全部位移坐标列向量x和第三步确定的控制力矩个数和位置，得到旋翼的振动控制微分方程；将旋翼振动控制微分方程改写成标准的状态空间方程表示形式；第三步中施加控制力矩位置处的节点位移作为输出响应向量y；

在有限元模型的基础上，基于状态空间的控制方法，实现对旋翼模型的振动主动控制。旋翼有限元模型的状态控制方程为：

$M_0\stackrel{··}{x}\left(t\right)+C_0\stackrel{·}{x}\left(t\right)+K_{0}x\left(t\right)=LU\left(t\right)+NF\left(t\right)---\left(1\right)$

其中，其中x＝[x₁,x₂,……,x_n]^T是有限元模型中的全部位移坐标列向量。有限元模型的自由度个数为n。这里M₀,C₀,K₀分别是旋翼有限元模型的质量矩阵，阻尼矩阵，刚度矩阵，均是n×n维。U(t)是p×1维的控制力向量，F(t)是m×1维的外部激励向量，这里的外部激励是气动载荷。L是n×p维的控制力定位矩阵，N是n×m维的外部激励定位矩阵。

把(1)式改写成：

$C_0\stackrel{·}{x}+M_0\stackrel{··}{x}+K_{0}x+0\stackrel{·}{x}=LU\left(t\right)+NF\left(t\right)---\left(2\right)$

同时有如下恒等式成立：

$M_0\stackrel{·}{x}+0\stackrel{··}{x}+0x+(-M_0)\stackrel{·}{x}=0---\left(3\right)$

可以联立写成：

$\left(\begin{array}{cc}{C}_0& M_0\\ M_0& 0\end{array}\right)\left(\begin{array}{c}\stackrel{·}{x}\\ \stackrel{··}{x}\end{array}\right)+\left(\begin{array}{cc}{K}_0& 0\\ 0& -M_0\end{array}\right)\left(\begin{array}{c}x\\ \stackrel{·}{x}\end{array}\right)=\left(\begin{array}{c}L\\ 0\end{array}\right)U\left(t\right)+\left(\begin{array}{c}N\\ 0\end{array}\right)F\left(t\right)---\left(4\right)$

通过逆矩阵的求解可知：

${\left(\begin{array}{cc}{C}_0& M_0\\ M_0& 0\end{array}\right)}^{-1}=\left(\begin{array}{cc}0& M_0^{-1}\\ M_0^{-1}& -M_0^{-1}{C}_0{M}_0^{-1}\end{array}\right)---\left(5\right)$

结合(4)式和(5)式可知：

$\left(\begin{array}{c}\stackrel{·}{x}\\ \stackrel{··}{x}\end{array}\right)=\left(\begin{array}{cc}0& I\\ -M_0^{-1}{K}_0& -M_0^{-1}{C}_0\end{array}\right)\left(\begin{array}{c}x\\ \stackrel{·}{x}\end{array}\right)+\left(\begin{array}{c}0\\ M_0^{-1}L\end{array}\right)U\left(t\right)+\left(\begin{array}{c}0\\ M^{-1}N\end{array}\right)F\left(t\right)---\left(6\right)$

由于压电纤维产生的是力矩作用，因而作用在旋翼上的控制力一定成对出现。假设共计作用有r个力矩，那么有：

p＝2r (7)

于是可以得到：

其中，u_i(t)与-u_i(t)成对出现组成力矩。

记并且记这里I是单位矩阵。

则(8)式可以表示成：

U(t)＝Tu(t) (9)

于是(6)式可以重新写成：

$\left(\begin{array}{c}\stackrel{·}{x}\\ \stackrel{··}{x}\end{array}\right)=\left(\begin{array}{cc}0& I\\ -M_0^{-1}{K}_0& -M_0^{-1}{C}_0\end{array}\right)\left(\begin{array}{c}x\\ \stackrel{·}{x}\end{array}\right)+\left(\begin{array}{c}0\\ M_0^{-1}LT\end{array}\right)u\left(t\right)+\left(\begin{array}{c}0\\ M_0^{-1}N\end{array}\right)F\left(t\right)---\left(10\right)$

其中的L是控制力列向量u(t)的定位矩阵。

$L=[l_1,l_2,...,l_r,l_{r+1},l_{r+2},...,l_{2r}]=\left(\begin{array}{cccccccc}{l}_{11}& l_{12}& ...& l_{1r}& l_{1,r+1}& l_{1,r+2}& ...& l_{1,2r}\\ l_{21}& l_{22}& ...& l_{2r}& l_{2,r+1}& l_{2,r+2}& ...& l_{2,2r}\\ .& .& & .& .& .& & .\\ .& .& & .& .& .& & .\\ .& .& & .& .& .& & .\\ l_{n1}& l_{n2}& ...& l_{nr}& l_{n,r+1}& l_{n,r+2}& ...& l_{n,2r}\end{array}\right)---\left(11\right)$

下面记z为2n×1维的状态列向量，

并且记：

$A=\left(\begin{array}{cc}0& I\\ -M_0^{-1}{K}_0& -M_0^{-1}{C}_0\end{array}\right)---\left(12\right)$

$B=\left(\begin{array}{c}0\\ M^{-1}LT\end{array}\right)---\left(13\right)$

$E=\left(\begin{array}{c}0\\ M_0^{-1}N\end{array}\right)---\left(14\right)$

于是(10)式可以写成：

$\stackrel{·}{z}\left(t\right)=Az\left(t\right)+Bu\left(t\right)+EF\left(t\right)---\left(15\right)$

设控制输出为位移列向量：

$y={[y_1,y_2,......,y_q]}^T=\left[\begin{array}{cc}H& 0\end{array}\right]·\left(\begin{array}{c}x\\ \stackrel{·}{x}\end{array}\right)=Cz---\left(16\right)$

其中H为控制输出的定位矩阵：

$H=\left(\begin{array}{c}{h}_1\\ h_2\\ .\\ .\\ .\\ h_q\end{array}\right)=\left(\begin{array}{cccc}{h}_{11}& h_{12}& ...& h_{1n}\\ h_{21}& h_{22}& ...& h_{2n}\\ .& .& & .\\ .& .& & .\\ .& .& & .\\ h_{q1}& h_{q2}& ...& h_{qn}\end{array}\right)---\left(17\right)$

且有：

C＝[H 0] (18)

把(16)式改写成：

y(t)＝Cz(t)+Du(t) (D＝0) (19)

联合(15)式和(19)式，考虑气动载荷作用下的旋翼振动控制问题可以表示成标准的状态空间描述如下。

$\left(\begin{array}{c}\stackrel{·}{z}\left(t\right)=Az\left(t\right)+Bu\left(t\right)+EF\left(t\right)\\ y\left(t\right)=Cz\left(t\right)+Du\left(t\right)\end{array}\right)---\left(20\right)$

A：系统矩阵；

B：控制输入矩阵；

C：输出矩阵；

D：直接传递矩阵，对一般的物理系统为0。

E：外部扰动定位矩阵；

这样就可以实现利用r个控制力实现对q个输出的控制。

第五步：根据第四步得到的状态空间方程，采用控制理论中的LQR算法，并且基于输出调节器的最优控制理论，进行旋翼振动的主动控制设计；

在输出调节器问题中，性能指标取为二次型形式如下：

$J=\frac{1}{2}{\int }_0^{\infty }[y{\left(t\right)}^{T}Qy\left(t\right)+u{\left(t\right)}^{T}Ru\left(t\right)]dt---\left(21\right)$

其中，利用第三步所有控制力矩中的力的大小组成列向量u和第四步输出响应向量y。最优控制理论要求J取最小值时的控制是最优控制。

控制理论中的LQR算法，要求矩阵Q和矩阵R是正定常数的权系数矩阵。计算中均取为单位对角矩阵。

最优控制要求J取最小值，即要求min J。u(t)不受约束，最优控制存在且唯一，由下式确定：

u(t)＝-R^-1B^TP·z(t)>

其中P是正定常数矩阵，且满足矩阵代数Riccati方程：

-PA-A^TP+PBR^-1B^TP-C^TQC＝0>

方程中矩阵P是未知的，其他矩阵均为已知。应用一般的商业数学软件，如MATLAB等，很容易通过求解矩阵代数Riccati方程得到未知矩阵P。根据矩阵Q和矩阵R以及第四步中的状态空间方程，解矩阵代数的Riccati方程，可以计算得到状态反馈增益矩阵G。

第六步：通过第五步中的LQR算法计算得到状态反馈增益矩阵G；通过G可以求出振动过程中施加的控制力矩的大小，从而开展旋翼振动主动控制；

根据矩阵Q和矩阵R以及第四步中的状态空间方程，解矩阵代数的Riccati方程，可以计算得到状态反馈增益矩阵G。

这时的状态反馈增益矩阵是：

G＝R^-1B^TP>

反馈控制向量为：

u＝-G·z(t) (25)

施加了反馈控制的状态运动轨迹是下列其次方程的解，即：

$\stackrel{·}{z}\left(t\right)=[A-BG]z\left(t\right)---\left(26\right)$

通过状态运动轨迹可以解出状态列向量z(t)。

这时的输出响应是：

y(t)＝Cz(t) (27)

第七步：基于第六步中旋翼实施振动主动控制之后的振动特性效果，判断进行主动控制之后的旋翼振动特性是否满足约束条件；如果满足，则结束；如果不满足，返回第三步，控制力矩的个数由原来的r变为r+1，然后重复第四步到第六步操作，直到旋翼的振动特性满足约束条件为止。

实施例：

针对某一型号旋翼，有限元模型见图2。旋翼长度2360mm，左端固支，右端为自由端。旋翼外表面铺设有压电纤维铺层，见图3，用以在通电情况下产生控制力矩。自由端受到F＝1000×sin(2.5×t)的外加激励，单位为N。

分别设置成无主动控制的振动和有LQR主动控制的振动，其中LQR控制中，权系数Q和R均设为单位阵。共计施加3个控制力矩，3个控制力矩的位置分别位于旋翼展向的25％、50％、75％处，示意图见图4。

首先进行无振动控制的振动响应分析，自由端位移响应历程见图5，图5中的实线表示旋翼端点关于时间的位移响应历程。图5中横坐标time代表时间历程，单位s；纵坐标displacement代表旋翼端点位移，单位mm。从图中可见最大位移约为6mm。

然后进行LQR振动主动控制。权系数Q和R均设为单位阵。自由端位移响应历程见图6。图6中的虚线表示旋翼端点关于时间的位移响应历程。图6中横坐标time代表时间历程，单位s；纵坐标displacement代表旋翼端点位移，单位mm。从图中可见最大位移约为1.5mm。

将有无振动控制的自由端位移响应曲线在图7中对比。图7中的实线表示无控制时的旋翼端点位移响应，虚线代表有控制时的旋翼端点位移响应。图7中横坐标time代表时间历程，单位s；纵坐标displacement代表旋翼端点位移，单位mm。可见采用LQR方法的振动主动控制可以有效降低旋翼的振动剧烈程度。

该实施例验证了本方法的有效性和适用性。

总之，本发明采用基于LQR算法的振动主动控制之后，旋翼的振动剧烈程度明显改善。本发明是主动控制理论与有限元分析的联合应用，便于在旋翼的设计阶段进行振动控制预估，提高了分析效率，通过采用主动控制的方法，改善了旋翼的振动特性，对于具体应用中有很强工程实践意义。

本发明未详细阐述部分属于本领域技术人员的公知技术。