CN103838141B - 一种面向控制的大型天线建模方法 - Google Patents

Abstract

本发明涉及一种面向控制的大型天线建模方法，其步骤是：第一步：ansvs模态分析；第二步：选择相关振型；第三步：计算模态参数；第四步：变换输入输出矩阵；第五步：分解各阶模态；第六步：计算各阶模态范数；第七步：模态缩聚；第八步：模态叠加。本发明的优点是：1)将天线柔性信息建立在控制模型中，使得模型更接近实际，在计算指向误差时也更为准确。2)刚性模型与柔性模型处于解耦状态，则模型分别输出刚性转角和柔性偏移，可分别对输出进行处理。

Description

一种面向控制的大型天线建模方法

技术领域

本发明属于天线技术领域，具体是一种面向控制的大型天线建模方法，用于实现将天线的柔性模态引入到伺服跟踪系统的模型中，提高分析精度。

背景技术

大型反射面天线以其照射效率高和边缘漏射低等优点，被广泛应用于卫星通讯、深空探索等领域，然而大型反射面天线指向控制器设计指标要求较高，我国“萤火工程”65米口径天线，其指向精度要求误差小于0.01°，只有精度达到要求，才能发挥天线口径效率，而随着口径的不断增大，天线结构会受到自重、风荷、雨雪等多种载荷作用，使得天线发生柔性变形，导致指向的偏差，而指向偏差又会引起天线电性能变差降低工作效率。当然随着天线工况的不同，这种影响也会随之变化。

对于如何降低外界载荷对天线指向控制的影响，除了需要先进控制理论的支持，还需要精确了解指向误差造成的因素，不光要克服由于传统的控制系统滞后对指向控制的阻碍，也需要解决由于天线自身结构柔性问题带来的指向误差，解决这一问题，就需要建立较为精确的面向控制的大型天线模型，然而传统的建模方法却忽略了影响天线指向的重要因素：天线自身柔性振动对指向的影响。

目前在国内外对提高天线指向的控制方法主要有以下几种：

1）通过不同的控制方法，包括PI控制、LQG控制及H无穷控制，解决由于各种外界环境所造成的指向误差，通过控制补偿使指向误差减小，通过比较三种控制方法，可知PI控制简单可靠但效果有限、而LQG控制和H无穷控制则大幅度提高了控制效果，同时结果表明通过控制方法来提高 指向精度作用有限，更需要在当前的天线驱动上加以改进。如在WodekGawronski.Antennacontrolsystems:FromPItoH∞,IEEEAntennasandPropagationMagazine,2001,vol.43,pp52-60.中就是用此种方法。

2）通过各种控制方法的组合，如速率环用PI控制、位置环用LQG控制，组合控制的指向精度为单一PI控制指向精度的二倍，同时若速率环为LQG控制、位置环为PID控制，则指向精度较单纯的PI控制提高十倍，若将LQG控制既应用于速率环又应用于位置环则效果是最好的，将会提高250倍。如在WodekGawronski.ControlSystemsofTheLargeMillimeterTelescope.IEEEAntennasandPropagationMagazine,2005,vol.47,pp41-49.中采用的就是这种方法。

3）主要通过分析了不同频域的风对天线转轴处及天线指向的影响，通过光学仪器测得在具体平均风速下，在不同风向和不同反射面位置处，该风造成的天线指向偏差，并建立图表，通过查表法对指向偏差进行平均补偿。如在NobuharuUkita,Wind-inducedpointingerrorsandsurface deformationofa10-msubmillimeterantenna,Conferenceon Ground-BasedandAirborneTelescopesШ,SanDiego,CA,2010中有所报导。

然而现有的提高指向精度的方法，要么将转轴处的刚性转角作为指向，分析外界载荷对其影响，并未考虑天线结构本身的柔性振动造成的指向偏差，降低了误差补偿准确率；要么利用光学仪器对不同工况下的指向误差数据建立表格，虽然这种方法相当于考虑了天线柔性振动，但其工作量大、效率低，且给予的补偿只是平均值，在一定程度上减小了指向误差，但效 果有限。

发明内容

本发明的目的在于避免上述现有技术的不足，提供一种面向控制的大型天线建模方法，将天线自身柔性信息考虑在建模过程中，克服传统建模对指向估计不准的缺点。

实现本发明目的的技术方案是：一种面向控制的大型天线建模方法，其步骤是：

第一步：以结构工况下反射面位置建立有限元模型，以某一频率F为约束条件，分析F以内的固有频率的振型，记为S1、S2......Sn；

第二步：以方位指向为模型的输出，第k个自由度为方位指向位移传感器位置，则振型其各阶模态对方位指向影响大小的初步估计取决于如下公式：

$C_{\mathrm{mq}}=C_{\mathrm{oq}}\mathrm{\Φ}=[00...1...0]\left[\begin{array}{cccc}{\mathrm{\φ}}_{{1n}_1}& {\mathrm{\φ}}_{{2n}_1}& \·\·\·& {\mathrm{\φ}}_{{\mathrm{nn}}_1}\\ {\mathrm{\φ}}_{{1n}_2}& {\mathrm{\φ}}_{{2n}_2}& \·\·\·& {\mathrm{\φ}}_{{\mathrm{nn}}_2}\\ \·\·\·& \·\·\·& \·\·\·& \·\·\·\\ {\mathrm{\φ}}_{{1n}_k}& {\mathrm{\φ}}_{{2n}_k}& \·\·\·& {\mathrm{\φ}}_{{\mathrm{nn}}_k}\\ \·\·\·& \·\·\·& \·\·\·& \·\·\·\\ {\mathrm{\φ}}_{{1n}_d}& {\mathrm{\φ}}_{{2n}_d}& \·\·\·& {\mathrm{\φ}}_{{\mathrm{nn}}_d}\end{array}\right]=[{\mathrm{\φ}}_{{1n}_k}{\mathrm{\φ}}_{{2n}_k}\·\·\·{\mathrm{\φ}}_{{\mathrm{nn}}_k}]$

我们称为第i阶模态对输出的增益值，

而柔性振荡引起的指向误差为：

e=C_mqq_m/l

其中q_m为柔性广义坐标，l为指向输出自由度到方位方向转轴的距离，若C_mq中最大元素为第k阶增益值虽然q_m未知，但当C_mq中的第i阶增益值远小于第k阶增益值根据经验有时，我们 认为第i阶模态对方位指向影响可忽略，同理俯仰方向准则一致；

第三步：提取所选各阶振型的参数，如各阶固有频率ω_i、各阶动能V_i及各阶第k个自由度振型向量φ_ik，根据建模原理计算柔性模型所需参数；

计算过程如下：

模态质量阵的提取，由于无法得到各自由度的质量而形成质量阵，因此理论上的模态质量无法通过下式得到：

M_m=Φ^TMΦ, (1)

以能量法简化计算各阶模态质量：

M_i=2V_i/(2πω_i)² (2)

其中，V_i为第i阶模态的总动能，ω_i为第i阶模态固有频率。

根据模态质量矩阵和固有频率矩阵计算模态刚度矩阵、模态阻尼矩阵和模态阻尼比矩阵：

Ω²=M_m ^-1K_m (3)

D_m=α₁K_m+α₂M_m (4)

$Z=0.5{M_m}^{-1}{D}_m{\mathrm{\Ω}}^{-1}=0.5{M_m}^{-\frac{1}{2}}{K_m}^{-\frac{1}{2}}{D}_m---\left(5\right)$

式中：Ω为自然频率矩阵，M_m称为模态质量阵，K_m称为模态刚度阵，D_m称为模态阻尼阵，Z为模态阻尼比矩阵。α₁,α₂为瑞丽阻尼系数,与结构固有频率相关:

α₁=(2(x₁ω₂-x₂ω₁)ω₁ω₂)/((ω₁+ω₂)(ω₂-ω₁))(6)

α₂=(2(x₂ω₂-x₁ω₁))/((ω₁+ω₂)(ω₂-ω₁))

其中ω₁,ω₂为结构前两阶固有频率，x₁,x₂一般取0.02。

第四步：再由模态振型矩阵和模态质量阵生成变换后的输入输出矩阵：

$\begin{array}{c}{B}_m={M_m}^{-1}{\mathrm{\Φ}}^T{B}_0,\\ C_{\mathrm{mq}}=C_{\mathrm{oq}}\mathrm{\Φ},\\ C_{\mathrm{mv}}=C_{\mathrm{ov}}\mathrm{\Φ}\end{array}---\left(7\right)$

式中Φ为柔性振型矩阵，B₀是输入矩阵，C_oq是位移输出矩阵，C_ov为速度输出矩阵。

第五步：根据模态建模方法各阶模态的可加性，将天线动力学方程分解改写为如下形式：

$\begin{array}{c}{\stackrel{\·\·}{q}}_{\mathrm{mi}}+2{\mathrm{\ξ}}_i{\mathrm{\ω}}_i{\stackrel{\·}{q}}_{\mathrm{mi}}+{{\mathrm{\ω}}_i}^2{q}_{\mathrm{mi}}=b_{\mathrm{mi}}u\\ y_i=c_{\mathrm{mqi}}{q}_{\text{mi}}{+c_{\mathrm{mvi}}\stackrel{\·}{q}}_{\mathrm{mi}},i=1,...,n,\\ y=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{y}_i\end{array}---\left(8\right)$

并将所得各阶模型改写为状态空间方程形式：

$A_{\mathrm{mi}}=\left[\begin{array}{cc}0& 1\\ -{{\mathrm{\ω}}_i}^2& -{2\mathrm{\ξ}}_i{\mathrm{\ω}}_i\end{array}\right],B_{\mathrm{mi}}=\left[\begin{array}{c}0\\ b_{\mathrm{mi}}\end{array}\right],C_{\mathrm{mi}}=[c_{\mathrm{mqi}}{c}_{\mathrm{mvi}}]---\left(9\right)$

第六步：根据缩聚原理，对初始状态空间方程中各阶模态进行范数计算：

$\left|\right|G_i|{|}_2\≅\frac{\left|\right|B_{\mathrm{mi}}\left|{|}_2\right|\left|C_{\mathrm{mi}}\right|{|}_2}{2\sqrt{{\mathrm{\ξ}}_i{\mathrm{\ω}}_i}}---\left(10\right)$

第七步：根据所得范数值将原状态空间方程进行分割，根据误差指标，截掉范数低的各阶模态，剩下的即为缩聚后的天线柔性模型：

误差定义为：

$e_2={\left(\underset{i=k+1}{\overset{n}{\mathrm{\Σ}}}{{\left|\right|G_i\left|\right|}_2}^2\right)}^{1/2}---\left(11\right)$

第八步：将所得天线柔性模型与刚性模型叠加生成面向控制的大型天线模型：

传统的模态法建模，对于刚性的处理忽略了结构外界阻尼D，将刚性模态简化为如下形式：

$\begin{array}{c}{\stackrel{\·\·}{q}}_{\mathrm{mi}}=b_{\mathrm{mi}}u\\ y_i=c_{\mathrm{mqi}}{q}_{\mathrm{mi}}+c_{\mathrm{mvi}}{\stackrel{\·}{q}}_{\mathrm{mi}}\end{array}---\left(12\right)$

从上式中可以看出，在恒定输入下，天线将以恒定加速度运动，这与实际是相违背的；因此，这里将引入外界阻尼D，天线刚性模型的转动惯量由ansys所得，外界阻尼D通过天线运动参数仿真估算得来，估算方法：在额定工作条件下，天线以某恒定的速度输出，调节阻尼D，使输出仿真值与实际值一致即可。

得到天线刚性模型：

$T=J\stackrel{\·\·}{\mathrm{\θ}}+D\stackrel{\·}{\mathrm{\θ}}---\left(13\right)$

柔性模型与刚性模型叠加，由于B₀与输入位置有关，对于不同的输入，比如控制力矩输入，和风力输入，B₀是不同的，因为两种输入激励的自由度位置不同，分别记为B₁和B₂。记控制力矩输入为T₁，风力输入为T₂。叠加过程如下：

${\mathrm{\Φ}}^T{B}_{1}u={\mathrm{\Φ}}^T{B}_1{T}_1={\mathrm{\Φ}}^T{B}_{1}J\stackrel{\·\·}{\mathrm{\θ}}+{\mathrm{\Φ}}^T{B}_{1}D\stackrel{\·}{\mathrm{\θ}}---\left(14\right)$

令：Φ^TB₁=B₁₁(15)

$B_{11}u=B_{11}{T}_1=B_{11}J\stackrel{\·\·}{\mathrm{\θ}}+B_{11}D\stackrel{\·}{\mathrm{\θ}}---\left(16\right)$

同理风力T₂，输入矩阵转换为B₁₂，最终得到天线模型为：

$\left[\begin{array}{cc}J& 0\\ B_{11}\·J& M_m\end{array}\right]\left[\begin{array}{c}\stackrel{\·\·}{\mathrm{\θ}}\\ {\stackrel{\·\·}{q}}_m\end{array}\right]+\left[\begin{array}{cc}D& 0\\ B_{11}\·D& D_m\end{array}\right]\left[\begin{array}{c}\stackrel{\·}{\mathrm{\θ}}\\ {\stackrel{\·}{q}}_m\end{array}\right]+\left[\begin{array}{cc}0& 0\\ 0& K_m\end{array}\right]\left[\begin{array}{c}\mathrm{\θ}\\ q_m\end{array}\right]=\left[\begin{array}{cc}{T}_1& 0\\ 0& T_2\end{array}\right]\left[\begin{array}{c}1\\ B_{12}\end{array}\right]---\left(17\right)$

本发明的优点是：

1)将天线柔性信息建立在控制模型中，使得模型更接近实际，在计算 指向误差时也更为准确。

2)刚性模型与柔性模型处于解耦状态，则模型分别输出刚性转角和柔性偏移，可分别对输出进行处理。

附图说明

下面结合实施例附图对本发明作进一步说明：

图1总流程图；

图2是第三步流程图分解；

图3振型三主反射面俯视振型图；

图4缩聚前后频域响应曲线；

图5刚性输出；

图6柔性输出。

具体实施方式

本发明的用到的模型有:1)有限元分析软件(ANSYS)中建立的天线有限元模型。

建模原理：

（1）初步参数提取原理：

首先根据已有的有限元建模方法，以结构工况下反射面位置建立有限元模型，以某一频率F为约束条件，一般情况下对结构的低频特性给予重视，因此F取100Hz，应用有限元软件分析F以内的各阶振型，假设共n阶，每一阶代号记为S1、S2......Sn，并选取n_d个自由度，n_d值一般大于2倍的阶数，根据振动理论，提取各阶振型，组成振型矩阵Φ，则Φ为n_d×n阶；

以方位指向为模型的输出，假设第k个自由度为方位指向位移传感器位置，则我们所关心的振型其各阶模态对方位指向影响大小的初步估计取决于如下公式：

$C_{\mathrm{mq}}=C_{\mathrm{oq}}\mathrm{\Φ}=[00...1...0]\left[\begin{array}{cccc}{\mathrm{\φ}}_{{1n}_1}& {\mathrm{\φ}}_{{2n}_1}& \·\·\·& {\mathrm{\φ}}_{{\mathrm{nn}}_1}\\ {\mathrm{\φ}}_{{1n}_2}& {\mathrm{\φ}}_{{2n}_2}& \·\·\·& {\mathrm{\φ}}_{{\mathrm{nn}}_2}\\ \·\·\·& \·\·\·& \·\·\·& \·\·\·\\ {\mathrm{\φ}}_{{1n}_k}& {\mathrm{\φ}}_{{2n}_k}& \·\·\·& {\mathrm{\φ}}_{{\mathrm{nn}}_k}\\ \·\·\·& \·\·\·& \·\·\·& \·\·\·\\ {\mathrm{\φ}}_{{1n}_d}& {\mathrm{\φ}}_{{2n}_d}& \·\·\·& {\mathrm{\φ}}_{{\mathrm{nn}}_d}\end{array}\right]=[{\mathrm{\φ}}_{{1n}_k}{\mathrm{\φ}}_{{2n}_k}\·\·\·{\mathrm{\φ}}_{{\mathrm{nn}}_k}]---\left(1\right)$

其中，C_oq为位移输出矩阵，阶数与自由度个数相同，哪个自由度放置位移传感器，则该自由度值为1，仿真中我们所关心的自由度输出即为1，其余为0，Φ为各阶振型矩阵，C_mq为模态输出矩阵。

而柔性振荡引起的指向误差为：

e=C_mqq_m/l (2)

其中q_m为柔性广义坐标，l为指向输出自由度到方位方向转轴的距离，若C_mq中最大元素为第k阶增益值虽然q_m未知，但当C_mq中的第i阶增益值远小于第k阶增益值根据经验有时，我们认为第i阶模态对方位指向影响可忽略，则此时我们将S1、S2......Sn中的第Si阶所有提取参数去掉，将剩下的各阶振型重新按固有频率由低到高排列记为第一阶、第二阶…第n阶，并重新改写振型矩阵Φ。同理俯仰方向准则一致；

（1）柔性建模原理：

结构的模态模型用模态坐标来表示，模态坐标各种参数由初步选取后的有限元软件中得来。此时我们如图2步骤301可得到各阶固有频率ω_i、 各阶动能V_i。此时自由度数仍为n_d，首先天线结构模型在节点坐标模型下可以被表达为如下二阶矩阵方程组：

$\begin{array}{c}M\stackrel{\·\·}{q}+D\stackrel{\·}{q}+\mathrm{Kq}=B_{0}u\\ y=C_{\mathrm{oq}}+C_{\mathrm{ov}}\stackrel{\·}{q}\·\end{array}---\left(3\right)$

在此方程中，q是n_d×1维的节点位移向量；是n_d×1维的节点速度向量；是n_d×1维的节点加速度向量；u是输入，y是输出；M是n_d×n_d维的质量矩阵；D是n_d×n_d维的阻尼矩阵；K是n_d×n_d维的刚度矩阵。输入矩阵B₀是n_d×1维，位移输出矩阵C_oq是1×n_d维，速度输出矩阵C_ov也是1×n_d维。质量矩阵是正定矩阵，刚度矩阵和阻尼矩阵均为半正定矩阵。

为了得到模态模型我们引入一个新的变量q_m，它满足如下关系：

q=Φq_m. (4)

Φ为柔性振型矩阵，q_m为柔性广义坐标，其中振型矩阵Φ为：

$\mathrm{\Φ}=[{\mathrm{\φ}}_1{\mathrm{\φ}}_2\·\·\·{\mathrm{\φ}}_n]=\left[\begin{array}{cccc}{\mathrm{\φ}}_{{1n}_1}& {\mathrm{\φ}}_{{2n}_1}& \·\·\·& {\mathrm{\φ}}_{{\mathrm{nn}}_1}\\ {\mathrm{\φ}}_{{1n}_2}& {\mathrm{\φ}}_{{2n}_2}& \·\·\·& {\mathrm{\φ}}_{{\mathrm{nn}}_2}\\ \·\·\·& \·\·\·& \·\·\·& \·\·\·\\ {\mathrm{\φ}}_{{1n}_k}& {\mathrm{\φ}}_{{2n}_k}& \·\·\·& {\mathrm{\φ}}_{{\mathrm{nn}}_k}\\ \·\·\·& \·\·\·& \·\·\·& \·\·\·\\ {\mathrm{\φ}}_{{1n}_d}& {\mathrm{\φ}}_{{2n}_d}& \·\·\·& {\mathrm{\φ}}_{{\mathrm{nn}}_d}\end{array}\right],---\left(5\right)$

把式(4)代入式(3)中，同时我们左乘Φ^T得到天线结构的柔性振动：

其中M_m称为模态质量阵，K_m称为模态刚度阵，D_m称为模态阻尼阵。

在实际应用中如图2步骤302，M_m=Φ^TMΦ,可根据能量法获得，具体如下式：M_i=2V_i/(2πω_i)² (7)

其中，V_i为第i阶模态的总动能，ω_i为第i阶模态固有频率。 计算图2步骤303各参数：D=α₁K+α₂M (8)

其中，α₁,α₂为瑞丽阻尼系数,与结构固有频率相关:

α₁=(2(x₁ω₂-x₂ω₁)ω₁ω₂)/((ω₁+ω₂)(ω₂-ω₁)) (9)

α₂=(2(x₂ω₂-x₁ω₁))/((ω₁+ω₂)(ω₂-ω₁))

其中ω₁,ω₂为结构前两阶固有频率，x₁,x₂一般取0.02。

再计算图2步骤304各参数：

K_m=Φ^TKΦ.D_m=Φ^TDΦ, (10)

将(6)式左乘M_m ^-1：

$\begin{array}{c}{\stackrel{\·\·}{q}}_m+2\mathrm{Z\Ω}{\stackrel{\·}{q}}_m+{\mathrm{\Ω}}^2{q}_m=B_{m}u,\\ y=C_{\mathrm{mq}}{q}_n+C_{\mathrm{mv}}{\stackrel{\·}{q}}_m\·\end{array}---\left(11\right)$

在上式中，Ω为自然频率矩阵，自然频率矩阵具有如下性质：

考虑到天线的自由振动没有外界阻尼情况下，

$M\stackrel{\·\·}{q}+\mathrm{Kq}=0.---\left(12\right)$

可得：

det(K-ω²M)=0.

Ω²=M_m ^-1K_m (13)

Z为模态阻尼比矩阵，也为对角线矩阵，

$\left[\begin{array}{cccc}{\mathrm{\ξ}}_1& 0& \·\·\·& 0\\ 0& {\mathrm{\ξ}}_2& \·\·\·& 0\\ \·\·\·& \·\·\·& \·\·\·& \·\·\·\\ 0& 0& \·\·\·& {\mathrm{\ξ}}_n\end{array}\right]---\left(14\right)$

其中ξ_i为第i阶模态阻尼比，Z矩阵我们可以通过M_m ^-1D_m=2ZΩ这个关系获得：

$Z=0.5{M_m}^{-1}{D}_m{\mathrm{\Ω}}^{-1}=0.5{M_m}^{-\frac{1}{2}}{K_m}^{-\frac{1}{2}}{D}_m---\left(15\right)$

同样模态输入B_m矩阵和模态输出矩阵C_mq,C_mv也存在这种变换：

$\begin{array}{c}{B}_m={M_m}^{-1}{\mathrm{\Φ}}^T{B}_0,\\ C_{\mathrm{mq}}=C_{\mathrm{oq}}\mathrm{\Φ},\\ C_{\mathrm{mv}}=C_{\mathrm{ov}}\mathrm{\Φ}\end{array}---\left(16\right)$

实际上，由于Ω和Z均为对角矩阵，因此式(11)也可以被写为如下形式：

$\begin{array}{c}{\stackrel{\·\·}{q}}_{\mathrm{mi}}+2{\mathrm{\ξ}}_i{\mathrm{\ω}}_i{\stackrel{\·}{q}}_{\mathrm{mi}}+{{\mathrm{\ω}}_i}^2{q}_{\mathrm{mi}}=b_{\mathrm{mi}}u\\ y_i=c_{\mathrm{mqi}}{q}_{\text{mi}}{+c_{\mathrm{mvi}}\stackrel{\·}{q}}_{\mathrm{mi}},i=1,...,n,\\ y=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{y}_i\end{array}---\left(17\right)$

其中，b_mi为B_m的第i行，c_mqi,c_mvi为C_mq,C_mv的第i列，y_i为系统的第i阶模态输出，需要注意的是系统的结构响应y为模态响应y_i的叠加，这正是使用模态坐标法得到系统性能的关键。

实际中，输入矩阵B₀与电机无关，它对应于结构的外载荷输入关系，其中中的元素除了在外载荷激励位置的自由度上为1，其他元素全为0，同样输出矩阵C_mq的元素除了在放置位移传感器的位置为1，其他元素也均为0。假如我们在有限元模型的第k个自由度处放置传感器，则有：

$C_{\mathrm{mq}}=C_{\mathrm{oq}}\mathrm{\Φ}=[00...1...0]\left[\begin{array}{cccc}{\mathrm{\φ}}_{{1n}_1}& {\mathrm{\φ}}_{{2n}_1}& \·\·\·& {\mathrm{\φ}}_{{\mathrm{nn}}_1}\\ {\mathrm{\φ}}_{{1n}_2}& {\mathrm{\φ}}_{{2n}_2}& \·\·\·& {\mathrm{\φ}}_{{\mathrm{nn}}_2}\\ \·\·\·& \·\·\·& \·\·\·& \·\·\·\\ {\mathrm{\φ}}_{{1n}_k}& {\mathrm{\φ}}_{{2n}_k}& \·\·\·& {\mathrm{\φ}}_{{\mathrm{nn}}_k}\\ \·\·\·& \·\·\·& \·\·\·& \·\·\·\\ {\mathrm{\φ}}_{{1n}_d}& {\mathrm{\φ}}_{{2n}_d}& \·\·\·& {\mathrm{\φ}}_{{\mathrm{nn}}_d}\end{array}\right]=[{\mathrm{\φ}}_{{1n}_k}{\mathrm{\φ}}_{{2n}_k}\·\·\·{\mathrm{\φ}}_{{\mathrm{nn}}_k}]---\left(18\right)$

这时，C_mq为传感器位置处各阶模态位移的组成，同样的对于在第k个自由度处放置电机激励，我们也得到输入矩阵的变化：

$B_m={M_m}^{-1}{\mathrm{\Φ}}^T{B}_0=\left[\begin{array}{c}{\mathrm{\φ}}_{{1n}_k}/m_{m1}\\ {\mathrm{\φ}}_{{2n}_k}/m_{m2}\\ \·\·\·\\ {\mathrm{\φ}}_{{\mathrm{nn}}_k}/m_{\mathrm{mm}}\end{array}\right]---\left(19\right)$

其中m_mi是第i阶模态质量，B_m为电机激励处模态位移与模态质量的组合。

同样我们也可以得到模型的状态空间表达式，从式（17）知，第i阶模态动力学方程可以被表示为如下形式的状态方程：

$\begin{array}{c}{\stackrel{\·}{q}}_{\mathrm{mi}}=A_i{x}_i+B_{i}u,\\ y_i=C_i{x}_i\\ y=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{y}_i\end{array}---\left(20\right)$

其中状态变量为： $x_i=\left\{\begin{array}{c}{q}_{\mathrm{mi}}\\ {\stackrel{\·}{q}}_{\mathrm{mi}}\end{array}\right\}---\left(21\right)$

相应的矩阵为：

$A_{\mathrm{mi}}=\left[\begin{array}{cc}0& 1\\ -{{\mathrm{\ω}}_i}^2& -{2\mathrm{\ξ}}_i{\mathrm{\ω}}_i\end{array}\right],B_{\mathrm{mi}}=\left[\begin{array}{c}0\\ b_{\mathrm{mi}}\end{array}\right],C_{\mathrm{mi}}=[c_{\mathrm{mqi}}{c}_{\mathrm{mvi}}]---\left(22\right)$

那么系统总模型的状态空间表达式就可以将各阶模态叠加起来，用

(A_m,B_m,C_m)表示：

$A=\mathrm{diag}\left(A_i\right)=\left[\begin{array}{cccccccc}\×& \×& 0& 0& \·\·\·& \·\·\·& 0& 0\\ \×& \×& 0& 0& \·\·\·& \·\·\·& 0& 0\\ 0& 0& \×& \×& \·\·\·& \·\·\·& 0& 0\\ 0& 0& \×& \×& \·\·\·& \·\·\·& 0& 0\\ \·\·\·& \·\·\·& \·\·\·& \·\·\·& \·\·\·& \·\·\·& \·\·\·& \·\·\·\\ \·\·\·& \·\·\·& \·\·\·& \·\·\·& \·\·\·& \·\·\·& \·\·\·& \·\·\·\\ 0& 0& 0& 0& \·\·\·& \·\·\·& \×& \×\\ 0& 0& 0& 0& \·\·\·& \·\·\·& \×& \×\end{array}\right],i=\mathrm{1,2},...,n,---\left(23\right)$

其中A_i为2×2阶的矩阵（非零元素用×表示），相应的有：

$B=\left[\begin{array}{c}{B}_1\\ B_2\\ \·\·\·\\ B_n\end{array}\right],C=[C_1{C}_2...C_n]---\left(24\right)$

（2）柔性模型缩聚原理：

通过模态截断缩聚：

在前面我们由模态坐标得到了系统状态空间方程模型，模型的状态变量可表达为如下形式：

$x=\left\{\begin{array}{c}{x}_1\\ x_2\\ \·\·\·\\ x_n\end{array}\right\}---\left(27\right)$

其中：

$x_i=\left\{\begin{array}{c}{q}_{\mathrm{mi}}\\ {\stackrel{\·}{q}}_{\mathrm{mi}}\end{array}\right\},---\left(28\right)$

我们令||G_i||₂表示第i阶模态的H₂范数，如下式表达：

$\left|\right|G_i|{|}_2\≅\frac{\left|\right|B_{\mathrm{mi}}\left|{|}_2\right|\left|C_{\mathrm{mi}}\right|{|}_2}{2\sqrt{{\mathrm{\ξ}}_i{\mathrm{\ω}}_i}}---\left(29\right)$

其中G_i我们定义为系统第i阶模态的传递函数：

$G_i\left(\mathrm{\ω}\right)=C_i{(\mathrm{j\ωI}-A_i)}^{-1}{B}_i=\frac{(c_{\mathrm{mqi}}+{\mathrm{j\ωc}}_{\mathrm{mvi}})b_{\mathrm{mi}}}{{{\mathrm{\ω}}_i}^2-{\mathrm{\ω}}^2+2j{\mathrm{\ξ}}_i{\mathrm{\ω}}_i\mathrm{\ω}},i=1,......,n,\left(30\right)$

模型各阶状态空间方程在前面已得到，如式（22）所示。B_mi和C_mi分别为第i阶模态的输入矩阵和输出矩阵，2范数定义为矩阵各元素平方和再开平方，那么||B_mi||₂和||C_mi||₂可以理解为输入增益和输出增益，ξ_i为第i阶模态阻尼比，ω_i为第i阶模态固有频率，那么||G_i||₂其实就是第i阶模态传递函数的增益，反映了在一定输入条件第i阶模态振幅的大小。

我们得到模型各阶模态的泛数后，将状态变量按照范数值从大到小从新排列，并作做如下分割：

$x=\left\{\begin{array}{c}{x}_r\\ x_t\end{array}\right\},---\left(31\right)$

其中x_r为我们要保留的状态向量，其范数值大，x_t为我们要截掉的状态 向量其范数值小，假设我们要保留的状态向量数为k，即x_r为2k阶，x_t为2(n-k)阶。我们也得到如下表达式：

$A_b=\left[\begin{array}{cc}{A}_r& 0\\ 0& A_t\end{array}\right],B_b=\left[\begin{array}{c}{B}_r\\ B_t\end{array}\right],C_b=\left[C_r{C}_t\right].---\left(32\right)$

降阶模型通过截掉矩阵A_b的最后2(n-k)行，2(n-k)列；截掉矩阵B_b的最后2(n-k)行；截掉矩阵C_b的最后2(n-k)列得到，即(A_r,B_r,C_r)。

截断标准由每阶模态的范数值衡量，那么究竟r与t又怎么取呢？这就要借助缩聚误差来判断。

H₂范数缩聚误差定义为：

e₂=||G-G_r||₂ (33)

其中G为未缩聚的传递函数模型，Gr为经过缩聚后的模型。前面提到过基于模态坐标的传递函数表达式可以表示为各个模态的叠加；因此，

$G=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{G}_i$ (34)

$G_r=\underset{i=1}{\overset{k}{\mathrm{\Σ}}}{G}_i$

设G-G_r=G_t，带入式（33），其中Gt为裁掉的部分。因此有，

e₂=||G_t||₂. (35)

又因为模态向量是平方可加的，我们便得到：

${\left|\right|G_t|{|}_2}^2=\underset{i=k+1}{\overset{n}{\mathrm{\Σ}}}{{\left|\right|G_i\left|\right|}_2}^2---\left(36\right)$

造成的误差为：

$e_2={\left(\underset{i=k+1}{\overset{n}{\mathrm{\Σ}}}{{\left|\right|G_i\left|\right|}_2}^2\right)}^{1/2}---\left(37\right)$

根据用途不同，我们缩聚准则也有区别，在某些情况下r已知，需要令 缩聚误差最小，那么只需要提取范数值最大的前r阶模态，同时可以计算缩聚模型的误差e₂，在另一些情况下e₂/||G||₂已知，则需要计算保留多少阶模态及哪些模态需要被保留比较合适，当然在误差允许的范围内，模型阶数越少越利于控制器的设计。

（3）多刚体建模：

由动力学方程所得：

$T=J\stackrel{\·\·}{\mathrm{\θ}}+D\stackrel{\·}{\mathrm{\θ}}---\left(38\right)$

其中J为天线惯量阵由ansys所得，D为外界阻尼阵。D由天线运动

参数仿真估算得来，方法为：在电机额定工作条件下，通过仿真调整阻

尼D使得天线输出（角速度）也为额定输出值。

（4）基于模态叠加的动力学建模：

为方便起见，将控制输入T₁用刚性转角表示，则有：

${\mathrm{\Φ}}^T{B}_{1}u={\mathrm{\Φ}}^T{B}_1{T}_1={\mathrm{\Φ}}^T{B}_{1}J\stackrel{\·\·}{\mathrm{\θ}}+{\mathrm{\Φ}}^T{B}_{1}D\stackrel{\·}{\mathrm{\θ}}---\left(39\right)$

在这里我们令：

Φ^TB₁=B₁₁

则有： $B_{11}u=B_{11}{T}_1=B_{11}J\stackrel{\·\·}{\mathrm{\θ}}+B_{11}D\stackrel{\·}{\mathrm{\θ}}---\left(40\right)$

前面已经介绍过B₀与输入位置有关，对于不同的输入，比如控制力矩输入，和风力输入，B₀是不同的，因为两种输入激励的自由度位置不同，我们记为B₁和B₂，同时模态输入矩阵也随之不同，我们记控制力矩模态输入矩阵为B₁₁，风力模态输入矩阵为B₁₂。

理论上控制力矩和风力输入都是均布在所接触得自由度上，但可做以简化，保证精度的同时，可以简化加载在几个主自由度上。

我们将所得柔性模型式（6）经过缩聚理论缩聚，截断不需要的各阶模态，将剩下的模型与刚性模型式（38）叠加，并将式（39）和式（40）代入，叠加后的模型：

$\left[\begin{array}{cc}J& 0\\ B_{11}\·J& M_m\end{array}\right]\left[\begin{array}{c}\stackrel{\·\·}{\mathrm{\θ}}\\ {\stackrel{\·\·}{q}}_m\end{array}\right]+\left[\begin{array}{cc}D& 0\\ B_{11}\·D& D_m\end{array}\right]\left[\begin{array}{c}\stackrel{\·}{\mathrm{\θ}}\\ {\stackrel{\·}{q}}_m\end{array}\right]+\left[\begin{array}{cc}0& 0\\ 0& K_m\end{array}\right]\left[\begin{array}{c}\mathrm{\θ}\\ q_m\end{array}\right]=\left[\begin{array}{cc}{T}_1& 0\\ 0& T_2\end{array}\right]\left[\begin{array}{c}1\\ B_{12}\end{array}\right]---\left(41\right)$

等式左边为输出，右边为输入，其中T1为电机输入力矩（均布力矩），T₂为风力（均布力）。

将本发明的建模方法应用于7.3米反射面天线进行仿真实验，建立有限元模型，反射面位置以俯仰角为45°为例，分析影响方位指向的柔性信息，分析100hz以内各阶振型，电机输入自由度简化为转轴加载力矩处2个对称自由度，反射面加载风力处5个自由度，以及以副反射面中心为输出指向处的自由度1个。

选取10阶初步建模，选取准则为对方位方向有指向影响的各阶振型，如图1所显示的振型主要发生在俯仰方向，影响俯仰指向比较明显，若以方位指向为模型的输出，则该阶模型的输出矩阵接近于[00]，意味着方位指向输出接近于0，所以此阶模态我们可以忽略。而图2显示的振型则主要发生在方位方向，因此我们在建立影响方位指向的模型时，该阶振型则需要被提取出来。观察图3所示振型可以发现，虽然主反射面并未发生明显振动，但副反射面及其支架振动明显，所以该阶模态虽然能量较小，但在初步建模过程中必须提取出来，当然该阶对方位方向指向影响有多大，后续缩聚模型计算范数后可得知。

对初步建立的10阶模型进行缩聚，以e₂/||G||₂<0.05为缩聚指标，得到保留模型应为6阶，且阶段误差为0.02，对比模型缩聚前后的频域响应曲线， 如图4所示，可以看出低频段模型基本保持一致。

通过ansys提取天线转动惯量J，同时通过simulink仿真，调试得出天线的外界阻尼D，调试参数为在额定电压110V时，该天线约以1°每秒的转速转动，并根据转动惯量J和阻尼D建立天线刚性模型，并与柔性模型叠加。

对所得天线模型进行仿真分析，将控制力矩和风扰加载与模型的输入端，观察其响应；将输入力矩先试做1000Nm，风扰力视为1000N的分析结果如图5、图6所示。

结果分析：从本发明的仿真实验中可以看出，由于柔性振动振幅相比刚性转角比较小，在传统控制中往往容易被忽略，所以极易导致指向精度不够准确，而随着大口径天线的发展，对指向的精度要求却越来越高，因此采用本发明的建模方法，可以在控制器设计阶段给予天线柔性信息以重视，提高天线指向精度。

本实施例没有详细叙述的部件和结构工艺属本行业的公知部件和常用结构或常用手段，这里不一一叙述。

Claims (1)

1.一种面向控制的大型天线建模方法，其步骤是：

第一步：以结构工况下反射面位置建立有限元模型，以某一频率F为约束条件，分析F以内的固有频率的振型，记为S1、S2......Sn；

第二步：以方位指向为模型的输出，第k个自由度为方位指向位移传感器位置，则振型其各阶模态对方位指向影响大小的初步估计取决于如下公式：

$C_{mq}=C_{oq}\mathrm{\&Phi;}=\left[\begin{array}{cccccc}0& 0& ...& 1& ...& 0\end{array}\right]\left[\begin{array}{cccc}{\mathrm{\&phi;}}_{1n_1}& {\mathrm{\&phi;}}_{2n_1}& ...& {\mathrm{\&phi;}}_{{\mathrm{nn}}_1}\\ {\mathrm{\&phi;}}_{1n_2}& {\mathrm{\&phi;}}_{2n_2}& ...& {\mathrm{\&phi;}}_{{\mathrm{nn}}_2}\\ ...& ...& ...& ...\\ {\mathrm{\&phi;}}_{1n_k}& {\mathrm{\&phi;}}_{2n_k}& ...& {\mathrm{\&phi;}}_{{\mathrm{nn}}_k}\\ ...& ...& ...& ...\\ {\mathrm{\&phi;}}_{1n_d}& {\mathrm{\&phi;}}_{2n_d}& ...& {\mathrm{\&phi;}}_{{\mathrm{nn}}_d}\end{array}\right]=\left[\begin{array}{cccc}{\mathrm{\&phi;}}_{1n_k}& {\mathrm{\&phi;}}_{2n_k}& ...& {\mathrm{\&phi;}}_{{\mathrm{nn}}_k}\end{array}\right]$

其中，称为第i阶模态对输出的增益值，C_mq为模态输出矩阵，而柔性振荡引起的指向误差为：

e＝C_mqq_m/l

其中q_m为柔性广义坐标，l为指向输出自由度到方位方向转轴的距离，若C_mq中最大元素为第k阶增益值虽然q_m未知，但当C_mq中的第i阶增益值远小于第k阶增益值根据经验有 时，认为第i阶模态对方位指向影响可忽略，同理俯仰方向准则一致；

第三步：提取所选各阶振型的参数，第i阶模态固有频率ω_i、第i阶模态的总动能V_i及各阶第k个自由度振型向量φ_ik，根据建模原理计算柔性模型所需参数；

计算过程如下：

模态质量阵的提取，由于无法得到各自由度的质量而形成质量阵，因此理论上的模态质量无法通过下式得到：

M_m＝Φ^TMΦ, (1)

以能量法简化计算各阶模态质量：

M_i＝2V_i/(2πω_i)² (2)

其中，V_i为第i阶模态的总动能，ω_i为第i阶模态固有频率；

根据模态质量矩阵和固有频率矩阵计算模态刚度矩阵、模态阻尼矩阵和模态阻尼比矩阵：

Ω²＝M_m ^-1K_m (3)

D_m＝α₁K_m+α₂M_m (4)

$Z=0.5{M_m}^{-1}{D}_m{\mathrm{\&Omega;}}^{-1}=0.5{M_m}^{-\frac{1}{2}}{K_m}^{-\frac{1}{2}}{D}_m---\left(5\right)$

式中：Ω为自然频率矩阵，M_m称为模态质量阵，K_m称为模态刚度阵，D_m称为模态阻尼阵，Z为模态阻尼比矩阵；α₁,α₂为瑞丽阻尼系数,与结构固有频率相关:

$\begin{array}{c}{\mathrm{\&alpha;}}_1=\left(2\left(x_1{\mathrm{\&omega;}}_2-x_2{\mathrm{\&omega;}}_1\right){\mathrm{\&omega;}}_1{\mathrm{\&omega;}}_2\right)/\left(\left({\mathrm{\&omega;}}_1+{\mathrm{\&omega;}}_2\right)\left({\mathrm{\&omega;}}_2-{\mathrm{\&omega;}}_1\right)\right)\\ {\mathrm{\&alpha;}}_2=\left(2\left(x_2{\mathrm{\&omega;}}_2-x_1{\mathrm{\&omega;}}_1\right)\right)/\left(\left({\mathrm{\&omega;}}_1+{\mathrm{\&omega;}}_2\right)\left({\mathrm{\&omega;}}_2-{\mathrm{\&omega;}}_1\right)\right)\end{array}---\left(6\right)$

其中ω₁,ω₂为结构前两阶固有频率，x₁,x₂一般取0.02；

第四步：再由模态振型矩阵和模态质量阵生成变换后的输入输出矩阵：

模态输入B_m矩阵和模态输出矩阵C_mq,C_mv存在这种变换：

$\begin{array}{c}{B}_m={M_m}^{-1}{\mathrm{\&Phi;}}^T{B}_0,\\ C_{mq}=C_{oq}\mathrm{\&Phi;},\\ C_{mv}=C_{ov}\mathrm{\&Phi;}\end{array}---\left(7\right)$

式中Φ为柔性振型矩阵，B₀是输入矩阵，C_oq是位移输出矩阵，C_ov为速度输出矩阵；

第五步：根据模态建模方法各阶模态的可加性，将天线动力学方程分解改写为如下形式：

$\begin{array}{c}{\stackrel{\&CenterDot;\&CenterDot;}{q}}_{mi}+2{\mathrm{\&xi;}}_i{\mathrm{\&omega;}}_i{\stackrel{\&CenterDot;}{q}}_{mi}+{{\mathrm{\&omega;}}_i}^2{q}_{mi}=b_{mi}u\\ y_i=c_{mqi}{q}_{mi}+c_{mvi}{\stackrel{\&CenterDot;}{q}}_{mi},i=1,...,n,\\ y=\underset{i=1}{\overset{n}{\mathrm{\&Sigma;}}}{y}_i\end{array}---\left(8\right)$

式中：b_mi为B_m的第i行，c_mqi,c_mvi为C_mq,C_mv的第i列，y_i为系统的第i阶模态输出，系统的结构响应y为模态响应y_i的叠加，ξ_i为第i阶模态阻尼比，

并将所得各阶模型改写为状态空间方程形式：

C_mi＝[c_mqi c_mvi](9)

第六步：根据缩聚原理，对初始状态空间方程中各阶模态进行范数计算：

$\left|\right|G_i|{|}_2\&cong;\frac{\left|\right|B_{mi}\left|{|}_2\right|\left|C_{mi}\right|{|}_2}{2\sqrt{{\mathrm{\&xi;}}_i{\mathrm{\&omega;}}_i}}---\left(10\right)$

第七步：根据所得范数值将原状态空间方程进行分割，根据误差指标，截掉范数低的各阶模态，剩下的即为缩聚后的天线柔性模型：

误差定义为：

$e_2={\left(\underset{i=k+1}{\overset{n}{\mathrm{\&Sigma;}}}\right|\left|G_i\right|{{|}_2}^2)}^{1/2}---\left(11\right)$

第八步：将所得天线柔性模型与刚性模型叠加生成面向控制的大型天线模型：

传统的模态法建模，对于刚性的处理忽略了结构外界阻尼D，将刚性模态简化为如下形式：

$\begin{array}{c}{\stackrel{\&CenterDot;\&CenterDot;}{q}}_{mi}=b_{mi}u\\ y_i=c_{mqi}{q}_{mi}+c_{mvi}{\stackrel{\&CenterDot;}{q}}_{mi}\end{array}---\left(12\right)$

从上式中可以看出，在恒定输入下，天线将以恒定加速度运动，这与实际是相违背的；因此，在这里将引入外界阻尼D，天线刚性模型的转动惯量由ansys所得，外界阻尼D通过天线运动参数仿真估算得来，估算方法：在额定工作条件下，天线以某恒定的速度输出，调节阻尼D，使输出仿真值与实际值一致即可；

得到天线刚性模型：

$T=J\stackrel{\&CenterDot;\&CenterDot;}{\mathrm{\&theta;}}+D\stackrel{\&CenterDot;}{\mathrm{\&theta;}}---\left(13\right)$

其中J为天线惯量阵，柔性模型与刚性模型叠加，由于B₀与输入位置有关，对于不同的控制力矩输入，和风力输入，B₀是不同的，因为两种输入激励的自由度位置不同，分别记为B₁和B₂；记控制力矩输入为T₁，风力输入为T₂；叠加过程如下：

${\mathrm{\&Phi;}}^T{B}_{1}u={\mathrm{\&Phi;}}^T{B}_1{T}_1={\mathrm{\&Phi;}}^T{B}_{1}J\stackrel{\&CenterDot;\&CenterDot;}{\mathrm{\&theta;}}+{\mathrm{\&Phi;}}^T{B}_{1}D\stackrel{\&CenterDot;}{\mathrm{\&theta;}}---\left(14\right)$

令：Φ^TB₁＝B₁₁(15)

$B_{11}u=B_{11}{T}_1=B_{11}J\stackrel{\&CenterDot;\&CenterDot;}{\mathrm{\&theta;}}+B_{11}D\stackrel{\&CenterDot;}{\mathrm{\&theta;}}---\left(16\right)$

同理风力T₂，输入矩阵转换为B₁₂，最终得到天线模型为：

$\left[\begin{array}{cc}J& 0\\ B_{11}\&CenterDot;J& M_m\end{array}\right]\left[\begin{array}{c}\stackrel{\&CenterDot;\&CenterDot;}{\mathrm{\&theta;}}\\ {\stackrel{\&CenterDot;\&CenterDot;}{q}}_m\end{array}\right]+\left[\begin{array}{cc}D& 0\\ B_{11}\&CenterDot;D& D_m\end{array}\right]\left[\begin{array}{c}\stackrel{\&CenterDot;}{\mathrm{\&theta;}}\\ {\stackrel{\&CenterDot;}{q}}_m\end{array}\right]+\left[\begin{array}{cc}0& 0\\ 0& K_m\end{array}\right]\left[\begin{array}{c}\mathrm{\&theta;}\\ q_m\end{array}\right]=\left[\begin{array}{cc}{T}_1& 0\\ 0& T_2\end{array}\right]\left[\begin{array}{c}1\\ B_{12}\end{array}\right]---\left(17\right).$