基于Lyapunov理论的卫星大挠性太阳能帆板分布式自适应振动抑制控制方法

Distributed adaptive vibration suppression control method of large solar panels for satellites based on Lyapunov theory

摘要

提出了一种分布式自适应振动抑制控制方法，有效抑制了卫星大挠性帆板在轨振动问题。考虑作动器和传感器的安装位置对控制效果的影响，通过分析邻接子模块之间振动和输出力的耦合影响，基于Lyapunov理论设计了完整的分布式自适应控制器。仿真结果表明，在外界持续干扰等情况下，本文设计的分布式控制系统与基于线性二次型调节器（LQR）的分布式控制方法相比抑振时间缩短了40%，与集中式控制方法相比抑振时间缩短了50%，得到了更良好的在轨抑振效果，对于卫星的稳定运行具有重要意义。

关键词

分布式控制 / 自适应控制 / 大挠性航天器 / 李雅普诺夫理论 / 振动抑制

Key words

distributed control / adaptive control / large flexible spacecraft / Lyapunov theory / vibration suppression

0 引言

太阳能帆板作为航天器的主要电力能源保障部件，目前正朝着大型化、挠性化方向发展^［1］。空间环境以及卫星姿态调整引起的耦合影响，使得帆板极易发生振动，且振动很难依靠被动抑振自行衰减，因此设计控制器进行主动振动控制对于卫星稳定运行具有重要意义。

针对太阳能帆板主动抑振控制的方法主要有两类：①集中式控制方法，中心计算机采集传感器信号并输出作动器控制力对应的电压信号；②分布式控制方法，由作动器传感器和控制电路组成控制节点，节点之间形成通信拓扑网络，实现分布式控制。Liu等^［2］利用拉格朗日公式推导了耦合动力学方程，设计了包含挠性附件在内的整星输出反馈控制器。Yuan等^［3］将改进的正位置反馈控制律应用于主动振动抑制，设计了基于耦合动力学的振动主动控制器。Luo等^［4］采用压电作动器实现了大型桁架结构的比例-微分（PD）振动控制。Wang等^［5］将时滞反馈控制器与PD控制器结合，在时域上研究了柔性航天器控制系统的稳定性，在频域上分析了柔性附件的阻尼作用。Rahman等^［6］采用比例-积分-微分（Proportion integration differentiation，PID）控制器，并基于智能层合板进行试验，证明该控制器在智能结构主动振动控制中的有效性。Tian等^［7］引入最优权重矩阵，设计了线性二次型调节器（Linear quadratic regulator，LQR）控制方案进行振动主动控制。Hu等^［8］设计了输出反馈控制器，能较好地抑制低频振动。苗双全等^［9］提出了一种输入成形与自适应滑模控制相结合的控制策略。上述控制方法均属于集中式控制方法，该类方法要求中心计算机的实时计算能力较强，但在控制节点失效时^［10］，系统容错性不高，而且太阳能帆板等部件尺寸成倍增加，集中式控制系统逐渐不能满足实际工程的需要。

Ji等^［11］针对一个中心体和两个对称挠性附件的航天器设计了分布式控制器，能够保证在执行机构局部故障时系统的稳定运行。Li等^［12］提出了基于分布式和分散观测器的容错控制方案，利用周期-小增益定理保证柔性航天器的闭环稳定性。Nakka等^［13］提出了一种用于航天器单自由度动态模型精确姿态控制的控制体系结构和算法，并设计了一种分布式驱动系统。Chen等^［14］基于分布式一致性协调控制理论，在无向或有向通信图下设计了无源控制器，可同时实现分布式姿态跟踪和振动抑制。上述分布式控制方法大多只考虑了邻接控制节点的振动耦合影响，有关邻接控制节点输出控制力的耦合影响因素需要进一步研究。

基于上述分析，针对分布式控制方法应用于航天器大挠性附件上的现状，考虑到执行机构输出控制力对于邻接子模块的影响，本文设计了一种分布式控制方法，并基于Lyapunov稳定性理论进行稳定性证明。运用Simulink进行仿真分析，与集中式控制方法及其他分布式控制算法进行对比分析，证明了本文方法的有效性，验证了控制律的正确性，并从实际角度出发，设置外界干扰情况，证明本文分布式控制方法具有更好的鲁棒性。

1 大挠性太阳能帆板子模块建模

以某卫星单侧太阳能帆板为研究对象，如图1所示。基于帆板建立坐标系，如图2所示。其中，大挠性帆板的三轴坐标系原点位于帆板与中心刚体链接的中心点处；x轴正方向为远离中心刚体方向；y轴与x轴垂直，且xy平面与帆板静止时的平面重合；z轴与x轴、y轴构成右手坐标系；根据项目需求中作动器/传感器对的数目，将大挠性太阳能帆板划分为2行4列的8个子模块，每个子模块中安装有一组作动器/传感器。

图2中， x _i （i=1，2，…，8）为帆板子模块振动的位移向量。本文忽略帆板在xoy平面的扭转，只考虑帆板在z轴方向上的振动。作动器安装在帆板上，根据控制需要输出控制力作用于帆板，使帆板实现主动振动抑制。作动器以两列的形式沿帆板坐标系的x轴方向排列，作动器采用压电陶瓷作动器，安装在大挠性太阳能帆板的表面。为避免运用压电元件时产生的控制溢出现象^［15］，传感器与作动器对位安装。帆板根据作动器的位置进行子模块划分。

压电陶瓷作动器/传感器是目前航天器挠性部件的主要抑振手段，针对如何利用有限资源提高振动抑制效果的问题，提出作动器/传感器位置优化配置方案，使得卫星大挠性帆板控制系统获得更好的控制效果。

分布式控制系统中的子模块耦合关系如图3所示。以子模块1为例，考虑子模块之间的耦合影响，需要注意的是，仅考虑与其边相邻接的两个子模块2和子模块i+1。其余子模块对其耦合影响视为作用于子模块2和子模块i+1之上，即其他子模块对该子模块的振动耦合影响和输出控制力耦合影响不做直接叠加计算。考虑到大挠性太阳能帆板的实际尺寸，该耦合影响较小，可以由其边界相邻的子模块表示。

本文分布式控制方法的整体设计思路是基于神经网络思想，将划分出的帆板子模块作为神经元；距离中心刚体最近的子模块1和子模块i+1既作为输入层又作为系统的输出层：作为输入层时，输入为激发函数，即所设计的控制律，作为输出层时，输出为子模块当前的状态信息；其余的子模块作为隐含层，且各个“神经元”之间的耦合关系项作为权重连接函数，整体形成类神经网络的系统结构。过程中放弃全局最优以实现子模块局部最优，使得距离中心刚体最近对其影响最大的子模块稳定，以保证挠性部件系统对航天器整体的振动影响最小。本文控制方案如图4所示。

对划分出的n个子模块进行动力学分析，为建立有限元模型，采用四节点板单元，利用哈密顿原理推导出板单元的动力学方程^［16］。通过组合可以得到控制子模块i的动力学方程，考虑前N阶模态，即有限元离散化的振动方程如下所示：

{\mathit{M}}_i {\ddot{\mathit{q}}}_i + {\mathit{C}}_i {\dot{\mathit{q}}}_i + {\mathit{K}}_i {\mathit{q}}_i = - {\mathit{B}}_{iu} {\mathit{u}}_i （1）

式中： {\mathit{q}}_i \in {\mathbf{R}}^{N\times \mathrm{1}}为子模块i的振动模态坐标； {\mathit{M}}_i \in {\mathbf{R}}^{N\times N}和 {\mathit{K}}_i \in {\mathbf{R}}^{N\times N}分别为子模块i的广义质量矩阵和刚度矩阵； {\mathit{C}}_i \in {\mathbf{R}}^{N\times N}为阻尼矩阵，且为瑞利阻尼，即 {\mathit{C}}_i = k_{\mathit{M}} {\mathit{M}}_{\mathit{i}} + k_{\mathit{K}} {\mathit{K}}_i， k_{\mathit{M}}、 k_{\mathit{K}}均为瑞利阻尼系数； {\mathit{B}}_{iu}为作动器的位置矩阵； {\mathit{u}}_i \in {\mathbf{R}}^{N\times \mathrm{1}}为控制输入，由压电陶瓷作动器加载电压后产生。

大挠性太阳能帆板整体的作动器方程和传感器方程的模态形式分别如下所示^［17］：

\begin{array}{c}{\ddot{\mathit{\eta }}}_{\mathrm{c}}(t)+\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left(\mathrm{2}{\xi }_i{\omega }_i\right){\dot{\mathit{\eta }}}_{\mathrm{c}}(t)+\mathrm{d}\mathrm{i}\mathrm{a}\mathrm{g}\left({\omega }_i^{\mathrm{2}}\right)·\\ {\mathit{\eta }}_{\mathrm{c}}(t)={\widehat{\mathit{B}}}_{\mathrm{c}}\mathit{u}(t)\end{array} （2）

\mathit{y} ( t ) = {\mathit{C}}_{\mathrm{a}} {\mathit{\Phi }}_{\mathrm{c}} {\ddot{\mathit{\eta }}}_{\mathrm{c}} ( t ) + {\mathit{C}}_{\mathrm{r}} {\mathit{\Phi }}_{\mathrm{c}} {\dot{\mathit{\eta }}}_{\mathrm{c}} ( t ) + {\mathit{C}}_{\mathrm{d}} {\mathit{\Phi }}_{\mathrm{c}} {\mathit{\eta }}_{\mathrm{c}} ( t ) （3）

式中： {\mathit{\eta }}_{\mathrm{c}} \in {\mathbf{R}}^{m\times \mathrm{1}}为前m阶主控模态； {\widehat{\mathit{B}}}_{\mathrm{c}}为 m \times n_c维的主控模态矩阵； {\omega }_i为第i阶固有频率； {\xi }_i为第i阶模态阻尼比； \mathit{y} \in {\mathbf{R}}^{n_o\times \mathrm{1}}为测量向量； {\mathit{C}}_{\mathrm{a}}、 {\mathit{C}}_{\mathrm{r}} \mathrm{、} {\mathit{C}}_{\mathrm{d}} \in {\mathbf{R}}^{n_{\mathrm{o}}\times n}分别为加速度、速度、位移传感器的输出矩阵； \mathit{u} ( t ) \in {\mathbf{R}}^{n_{\mathrm{c}}\times \mathrm{1}}为作动器所产生的控制力向量； {\mathit{\Phi }}_{\mathrm{c}}为m阶主控模态振型组成的 n \times m维矩阵； {\mathit{\varphi }}_i = {\left[{\varphi }_{\mathrm{1}i},{\varphi }_{\mathrm{2}i},\cdots ,{\varphi }_{ni}\right]}^{\mathrm{T}}为第i阶模态振型。

2 作动器传感器位置优化配置

为了得到更好的控制效果，需要在各个子模块中安装一对作动器/传感器进行位置优化配置。考虑到航天器在轨的复杂工作环境，为了避免振动初值及控制律的影响，运用基于系统可控/可观度的优化准则方法进行寻优计算。

根据上述太阳能帆板的模态形式方程（2）（3），可控可观度的描述需要列写系统的状态方程和输出方程，所以取 \mathit{x} = {\left({\dot{\eta }}_{\mathrm{1}},{\omega }_{\mathrm{1}}{\eta }_{\mathrm{1}},\cdots ,{\dot{\eta }}_n,{\omega }_n{\eta }_n\right)}^{\mathrm{T}}，得状态方程为：

\dot{\mathit{x}} ( t ) = \tilde{\mathit{A}} \mathit{x} ( t ) + \tilde{\mathit{B}} \mathit{u} ( t ) （4）

式中： \tilde{\mathit{A}} = \mathrm{b} \mathrm{l} \mathrm{k} \mathrm{d} \mathrm{i} \mathrm{a} \mathrm{g} ( {\mathit{A}}_i )； \tilde{\mathit{B}} = ( {\mathit{B}}_{\mathrm{1}}^{\mathrm{T}} , {\mathit{B}}_{\mathrm{2}}^{\mathrm{T}} , \cdots , {\mathit{B}}_n^{\mathrm{T}} {)}^{\mathrm{T}}； {\mathit{A}}_i = \left(\begin{array}{cc}-\mathrm{2}{\zeta }_i{\omega }_i& -{\omega }_i\\ {\omega }_i& \mathrm{0}\end{array}\right)； {\mathit{B}}_{\mathrm{i}} = \left(\begin{array}{l}{\mathit{b}}_i^{\mathrm{T}}\\ \mathrm{0}\end{array}\right)， {\mathit{b}}_i^{\mathrm{T}}为 {\widehat{\mathit{B}}}_{\mathrm{c}}的第i行。

输出方程为：

\mathit{y} ( t ) = \tilde{\mathit{C}} \mathit{x} ( t ) + \tilde{\mathit{D}} \mathit{u} ( t ) （5）

式中： \tilde{\mathit{C}} = ( {\mathit{C}}_{\mathrm{1}} , {\mathit{C}}_{\mathrm{2}} , \cdots , {\mathit{C}}_n )和 \tilde{\mathit{D}}取值为 {\mathit{C}}_i = \left(\begin{array}{cc}{c}_{r_i}& \mathrm{0}\end{array}\right)， \tilde{\mathit{D}} = \mathrm{0}；传感器选择速度传感器。

结构特性决定系统稳定，当终止时间 T \to \mathrm{\infty }，可控性Gram矩阵 {\mathit{W}}_{\mathrm{c}} ( \mathrm{0} , T ) \to {\mathit{W}}_{\mathrm{c}}，可观性Gram矩阵 {\mathit{W}}_{\mathrm{o}} ( \mathrm{0} , T ) \to {\mathit{W}}_{\mathrm{o}}，分别满足如下Lyapunov方程：

\tilde{\mathit{A}} {\mathit{W}}_{\mathrm{c}} + {\mathit{W}}_{\mathrm{c}} {\tilde{\mathit{A}}}^{\mathrm{T}} + \tilde{\mathit{B}} {\tilde{\mathit{B}}}^{\mathrm{T}} = \mathrm{0} （6）

{\tilde{\mathit{A}}}^{\mathrm{T}} {\mathit{W}}_{\mathrm{o}} + {\mathit{W}}_{\mathrm{o}} \tilde{\mathit{A}} + {\tilde{\mathit{C}}}^{\mathrm{T}} \tilde{\mathit{C}} = \mathrm{0} （7）

从式（6）中可以看出：可控性Gram矩阵 {\mathit{W}}_{\mathrm{c}}相关的矩阵方程（6）中 \tilde{\mathit{B}}矩阵包含了作动器的位置信息，该Lyapunov方程求解得到的可控性Gram矩阵也同样与作动器的位置密切相关；同理，可观性Gram矩阵 {\mathit{W}}_{\mathrm{o}}也与传感器的位置密切相关。求解式（6）可以得到可控性Gram矩阵 {\mathit{W}}_{\mathrm{c}}，利用 \tilde{\mathit{A}}的块对角形式，将 {\mathit{W}}_{\mathrm{c}}拆分成为两行、两列的块矩阵 {\mathit{W}}_{\mathrm{c}ij}，故式（6）可以变换为 {\mathit{A}}_i {\mathit{W}}_{\mathrm{c}ij} + {\mathit{W}}_{\mathrm{c}ij} {\mathit{A}}_j^{\mathrm{T}} + {\mathit{B}}_i {\mathit{B}}_j^{\mathrm{T}} = \mathrm{0}。得到可控性Gram矩阵的解析形式为：

{\mathit{W}}_{\mathrm{c}ij} = \frac{{\mathit{\beta }}_{ij}}{{\mathit{d}}_{ij}} · \left(\begin{array}{cc}\mathrm{2}{\omega }_i{\omega }_j\left({\zeta }_j{\omega }_i+{\zeta }_i{\omega }_j\right)& {\omega }_j\left({\omega }_j^{\mathrm{2}}-{\omega }_i^{\mathrm{2}}\right)\\ -{\omega }_i\left({\omega }_j^{\mathrm{2}}-{\omega }_i^{\mathrm{2}}\right)& \mathrm{2}{\omega }_i{\omega }_j\left({\zeta }_i{\omega }_i+{\zeta }_j{\omega }_j\right)\end{array}\right) （8）

式中： {\mathit{\beta }}_{ij} = {\mathit{b}}_i^{\mathrm{T}} {\mathit{b}}_j，其中， {\mathit{b}}_i^{\mathrm{T}}为取 \mathit{x} = ( {\dot{\eta }}_{\mathrm{1}} , {\omega }_{\mathrm{1}} {\eta }_{\mathrm{1}} , \cdots , {\dot{\eta }}_n , {\omega }_n {\eta }_n {)}^{\mathrm{T}}作为作动器传感器方程状态变量时状态方程 B 矩阵的第i行，n为帆板离散化后的自由度； d_{ij} = \mathrm{4} {\omega }_i {\omega }_j \left({\zeta }_i{\omega }_i+{\zeta }_j{\omega }_j\right) \left({\zeta }_j{\omega }_i+{\zeta }_i{\omega }_j\right) + {\left({\omega }_j^{\mathrm{2}}-{\omega }_i^{\mathrm{2}}\right)}^{\mathrm{2}}。

此时的Gram矩阵解析形式相比于直接求解Lyapunov方程计算量极大减少。

同理，求解式（7）可以得到可观性Gram矩阵 {\mathit{W}}_{\mathrm{o}}，将 {\mathit{W}}_{\mathrm{o}}矩阵拆分成两行、两列的块矩阵，则计算后可得到 {\mathit{W}}_{\mathrm{o}}的解析形式为：

{\mathit{W}}_{\mathrm{o}ij} = \frac{{\mathit{\gamma }}_{ij}}{{\mathit{d}}_{ij}} · \left(\begin{array}{cc}\mathrm{2}{\omega }_i{\omega }_j\left({\zeta }_j{\omega }_i+{\zeta }_i{\omega }_j\right)& -{\omega }_j\left({\omega }_j^{\mathrm{2}}-{\omega }_i^{\mathrm{2}}\right)\\ {\omega }_i\left({\omega }_j^{\mathrm{2}}-{\omega }_i^{\mathrm{2}}\right)& \mathrm{2}{\omega }_i{\omega }_j\left({\zeta }_i{\omega }_i+{\zeta }_j{\omega }_j\right)\end{array}\right) （9）

式中： {\mathit{\gamma }}_{ij} = {\mathit{c}}_{ri}^{\mathrm{T}} {\mathit{c}}_{rj}， {\mathit{c}}_{ri}^{}为传感器为速度传感器的输出方程 C 矩阵的第i列；观测量为速度输出。

根据上述计算，已知系统的可控性和可观性Gram矩阵的解析形式，可以计算得到Gram矩阵的奇异值，奇异值代表了主控模态的可控度和可观度。对于帆板这种弱阻尼密集模态系统，提高或降低主控模态中一些模态的可控度或可观度会引起其他模态可控度可观度相应的升高或下降。

本文采用的作动器传感器优化配置准则为：

J_{\mathrm{c}} = \left(\frac{\mathrm{1}}{\mathrm{2}\times m}{\displaystyle \sum _{i=\mathrm{1}}^{\mathrm{2}\times m}}{\mathit{\sigma }}_{\mathrm{c}i}\right) \sqrt[\mathrm{2}\times m]{{\displaystyle \prod _{i=\mathrm{1}}^{\mathrm{2}\times m}}{\mathit{\sigma }}_{\mathrm{c}i}} （10）

式中：m为主控模态数； {\mathit{\sigma }}_{\mathrm{c}i}为可控性Gram矩阵 {\mathit{W}}_{\mathrm{c}}的奇异值，该值代表所有主控模态的可控度。

优化准则乘积分为两部分：前一部分为所有奇异值的算术平均值，该项的意义在于使得所有奇异值最大；后一部分是所有奇异值的几何平均值，该项的意义在于使得所有的奇异值集中分布尽可能最大，避免出现有某些奇异值很大导致平均值虽然最大，但是奇异值分布不集中、有较小奇异值出现的情况。因此，作动器的优化配置准则最大时，可以保证整个卫星太阳能帆板系统的振动抑制效果最佳。同理，本文采用的传感器优化配置准则为：

J_{\mathrm{o}} = \left(\frac{\mathrm{1}}{\mathrm{2}\times m}{\displaystyle \sum _{i=\mathrm{1}}^{\mathrm{2}\times m}}{\mathit{\sigma }}_{\mathrm{o}i}\right) \sqrt[\mathrm{2}\times m]{{\displaystyle \prod _{i=\mathrm{1}}^{\mathrm{2}\times m}}{\mathit{\sigma }}_{\mathrm{o}i}} （11）

式中： {\mathit{\sigma }}_{\mathrm{o}i}为可观性Gram矩阵 {\mathit{W}}_{\mathrm{o}}的奇异值。

与式（10）相同，该传感器的位置优化配置准则可以使得所有主控模态的可观度集中分布且达到最大。

在作动器/传感器位置优化过程中，可以将式（8）和式（9）的可控性Gram矩阵和可观性Gram矩阵代入式（10）和式（11）优化准则中，这样可以减少计算量，避免求解高阶Lyapunov方程。

帆板的各项参数如表1所示。经过上文的式（10）（11）可以计算出相应的可控度准则系数J _c和可观度准则系数J _o，其结果数量级相差巨大，要得出系统最终的可控可观度综合优化准则J，需要对原有两种优化准则结果进行归一化处理。将两个优化准则J _c 、J _o转化为同一数量级，即将J _c乘以1000化为10²数量级，同理，将J _o除以10，得到同一数量级的 J_{\mathrm{c}}^{\mathrm{\text{'}}} \mathrm{、} J_{\mathrm{o}}^{\mathrm{\text{'}}}，即 J_{\mathrm{c}}^{\mathrm{\text{'}}} = \mathrm{1000} \times J_{\mathrm{c}}^{}， J_{\mathrm{o}}^{\mathrm{\text{'}}} = \mathrm{0.1} \times J_{\mathrm{o}}^{}。

由于作动器和传感器是对位安装，故把两个相同数量级的准则系数合成可控可观度综合优化准则J： J = \mathrm{1} / ( {\alpha }_{\mathrm{1}} J_{\mathrm{c}}^{\mathrm{\text{'}}} + {\alpha }_{\mathrm{2}} J_{\mathrm{o}}^{\mathrm{\text{'}}} )。其中， {\alpha }_{\mathrm{1}} + {\alpha }_{\mathrm{2}} = \mathrm{1}， {\alpha }_{\mathrm{1}}为可控度的权重， {\alpha }_{\mathrm{2}}为可观度的权重，本文中 {\alpha }_{\mathrm{1}} = {\alpha }_{\mathrm{2}} = \mathrm{0.5}，即 J = \mathrm{1} / ( \mathrm{0.5} J_{\mathrm{c}}^{\mathrm{\text{'}}} + \mathrm{0.5} J_{\mathrm{o}}^{\mathrm{\text{'}}} )。

运用遗传算法进行目标函数的寻优计算，可以得到如表2所示的8个坐标点结果，在太阳能帆板上的位置如图5所示。

3 考虑控制力耦合的分布式控制方法设计

基于大挠性太阳能帆板子模块的振动方程设计分布式控制律，需要考虑两方面的耦合影响因素。如图6所示，帆板受到外界空间环境的干扰，产生难以自行衰减的振动。对其进行主动振动控制，以相邻接的序号为i和j的两个子模块为研究对象，说明邻接太阳能帆板子模块之间的耦合影响。

（1）邻接子模块j振动产生的干扰力矩对本子模块i的耦合影响。与航天器刚柔耦合情形中的挠性附件振动与中心刚体平动的耦合关系类比。

设定此一项的协同控制项为耦合系数 {\mathit{C}}_{i{n}_k}（其中， i为控制子模块序号， n_k为第k个邻接子模块序号）乘以邻接子模块模态坐标的二阶导数 {\ddot{\mathit{q}}}_{{\mathit{n}}_k}，即 {\mathit{C}}_{i{n}_k} {\ddot{\mathit{q}}}_{{\mathit{n}}_k}。

（2）邻接子模块产生的控制力对本子模块的耦合影响。当子模块根据所设计的控制律产生相应控制力时，该力矩不仅会在当前时刻t作用到自身的振动抑制上，还会在下一时刻t+1作用在邻接子模块j上。所以，子模块j需要对自身的控制律添加一个协同项，以保证在t+1时刻子模块所产生的控制力能够抵消邻接子模块的控制力影响，所以需要设定此一项的协同控制项为控制力耦合系数乘以邻接子模块的控制力，即 k_{n_k} {\mathit{f}}_{a{n}_k}。

状态反馈控制力由作动器产生，即 {\mathit{f}}_{ai} = {\mathit{B}}_{iu} {\mathit{u}}_i， {\mathit{f}}_{ai} = [ f_{i\mathrm{1}} ( t ) , f_{i\mathrm{2}} ( t ) , \cdots , f_{iN} {(t)]}^{\mathrm{T}}。设置分布式控制律为：

\begin{array}{l}{\mathit{f}}_{ij}(t)=-\left[\alpha {\mathit{q}}_{ij}\right.+\beta {\dot{\mathit{q}}}_{ij}+{\gamma }_{n_k}\sum \left({\dot{\mathit{q}}}_{n_{k}j}\right.{\ddot{\mathit{q}}}_{ij}^T{\dot{\mathit{q}}}_{n_{k}j}+\\ {\dot{\mathit{q}}}_{n_{k}j}\left.{\dot{\mathit{q}}}_{ij}^{\mathrm{T}}{\ddot{\mathit{q}}}_{n_{k}j}\right)+k_{n_k}\sum \left({\mathit{f}}_{a{n}_k}\right.{\ddot{\mathit{q}}}_{n_{k}j}^T{\mathit{f}}_{a{n}_k}+\\ {\mathit{f}}_{a{n}_k}{\ddot{\mathit{q}}}_{n_{k}j}^{\mathrm{T}}\left.\left.{\dot{\mathit{f}}}_{a{n}_k}\right)\right]\end{array} （12）

式中： \alpha = - R \le \mathrm{0}，R为任意正常数； ； {\gamma }_{n_k} = - {\mathit{C}}_{i{n}_{k}j}；前两项 \alpha {\mathit{q}}_{ij} + \beta {\dot{\mathit{q}}}_{ij}为状态反馈控制器； {\gamma }_{n_k} \sum ( {\dot{\mathit{q}}}_{n_{k}j} {\ddot{\mathit{q}}}_{ij}^{\mathrm{T}} {\dot{\mathit{q}}}_{n_{k}j} + {\dot{\mathit{q}}}_{n_{k}j} {\dot{\mathit{q}}}_{ij}^{\mathrm{T}} {\ddot{\mathit{q}}}_{n_{k}j}）为邻接子模块的振动耦合影响项； k_{n_k} \sum ( {\mathit{f}}_{a{n}_k} {\ddot{\mathit{q}}}_{n_{k}j}^{\mathrm{T}} {\mathit{f}}_{a{n}_k} + {\mathit{f}}_{a{n}_k} {\ddot{\mathit{q}}}_{n_{k}j}^{\mathrm{T}} {\dot{\mathit{f}}}_{a{n}_k} )为邻接子模块的作动器输出控制力耦合影响项。

与其他控制律相比，所设计的控制律考虑了邻接子模块作动器输出控制力的耦合影响，使得控制过程中的控制力计算更为精确，减少抑振过程中的振幅。所以，基于分布式控制方法的太阳能帆板子模块振动抑制方程可以表示为：

\begin{array}{l}{M}_{ij}{\ddot{\mathit{q}}}_{ij}+C_{ij}{\dot{\mathit{q}}}_{ij}+K_{ij}{\mathit{q}}_{ij}=\alpha {\mathit{q}}_{ij}+\beta {\dot{\mathit{q}}}_{ij}+\\ {\gamma }_{n_k}\sum \left({\dot{\mathit{q}}}_{n_{k}j}{\ddot{\mathit{q}}}_{ij}^T{\dot{\mathit{q}}}_{n_{k}j}+{\dot{\mathit{q}}}_{n_{k}j}{\dot{\mathit{q}}}_{ij}^T{\ddot{\mathit{q}}}_{n_{k}j}\right)+\\ k_{n_k}\sum \left({\mathit{f}}_{a{n}_k}{\ddot{\mathit{q}}}_{n_{k}j}^T{\mathit{f}}_{a{n}_k}+{\mathit{f}}_{a{n}_k}{\ddot{\mathit{q}}}_{n_{k}j}^T{\dot{\mathit{f}}}_{a{n}_k}\right)\end{array} （13）

为验证整个系统的稳定性，首先引出以下关于Lyapunov稳定性理论^［18］。

定理1 假设系统的状态方程： \dot{\mathit{x}} = \mathit{f} ( \mathit{x} , t )， \mathit{f} ( \mathrm{0} , t ) = \mathrm{0} , \forall t。如果存在一个具有连续偏导数的标量函数 V ( \mathit{x} , t )，并且满足条件：① V ( \mathit{x} , t )是正定的；② \dot{V} ( \mathit{x} , t )是负定的。那么系统是渐进稳定的。

然后，选取如下的Lyapunov函数：

\begin{array}{l}V=\left[{\dot{\mathit{q}}}_{ij}^T{M}_{ij}{\dot{\mathit{q}}}_{ij}\right.+{\mathit{q}}_{ij}^{\mathrm{T}}{K}_{ij}{\mathit{q}}_{ij}+{\mathit{q}}_{ij}^{\mathrm{T}}R{\mathit{q}}_{ij}+\\ \sum C_{i{n}_{k}j}({\dot{\mathit{q}}}_{ij}^T{\dot{\mathit{q}}}_{n_{k}j}{)}^{\mathrm{2}}+\\ \sum k_{n_k}\left.({\dot{\mathit{q}}}_{ij}^T{\mathit{f}}_{a{n}_k}{)}^{\mathrm{2}}\right]/\mathrm{2}\end{array} （14）

显然，式（14）恒大于0，且对于时间可求导。 式（14）对时间进行一阶求导可得：

\begin{array}{l}{V}^{\mathrm{\text{'}}}={\dot{\mathit{q}}}_{ij}^{\mathrm{T}}{M}_{ij}{\ddot{\mathit{q}}}_{ij}+{\mathit{q}}_{ij}^{\mathrm{T}}{K}_{ij}{\dot{\mathit{q}}}_{ij}+{\mathit{q}}_{ij}^{\mathrm{T}}R{\dot{\mathit{q}}}_{ij}+\\ \sum C_{i{n}_{k}j}({\dot{\mathit{q}}}_{ij}^T{\dot{\mathit{q}}}_{n_{k}j})({\ddot{\mathit{q}}}_{ij}^{\mathrm{T}}{\dot{\mathit{q}}}_{n_{k}j}+{\dot{\mathit{q}}}_{ij}^{\mathrm{T}}{\ddot{\mathit{q}}}_{n_{k}j})+\\ \sum k_{n_k}({\dot{\mathit{q}}}_{ij}^T{\mathit{f}}_{a{n}_k})({\ddot{\mathit{q}}}_{ij}^{\mathrm{T}}{\mathit{f}}_{a{n}_k}+{\dot{\mathit{q}}}_{ij}^{\mathrm{T}}{\dot{\mathit{f}}}_{a{n}_k})\end{array} （15）

式中： {\dot{\mathit{q}}}_{ij}项由上述的子模块振动方程得出。

\begin{array}{l}{\ddot{\mathit{q}}}_{ij}=-M_{ij}^{-\mathrm{1}}{C}_{ij}{\dot{\mathit{q}}}_{ij}-M_{ij}^{-\mathrm{1}}{K}_{ij}{\mathit{q}}_{ij}-\\ M_{ij}^{-\mathrm{1}}{B}_{iuj}{\mathit{u}}_{ij}\end{array} （16）

代入Lyapunov函数的一阶导数可得：

\begin{array}{l}{\mathit{V}}^{\mathrm{\text{'}}}={\dot{\mathit{q}}}_{ij}^{\mathrm{T}}{M}_{ij}(-M_{ij}^{-\mathrm{1}}{C}_{ij}{\dot{\mathit{q}}}_{ij}-M_{ij}^{-\mathrm{1}}{K}_{ij}{\mathit{q}}_{ij}-\\ M_{ij}^{-\mathrm{1}}{B}_{iuj}{\mathit{u}}_{ij})+{\mathit{q}}_{ij}^{\mathrm{T}}{K}_{ij}{\dot{q}}_{ij}+{\mathit{q}}_{ij}^{\mathrm{T}}R{\dot{\mathit{q}}}_{ij}+\\ \sum C_{i{n}_{k}j}({\dot{\mathit{q}}}_{ij}^T{\dot{\mathit{q}}}_{n_{k}j})({\ddot{\mathit{q}}}_{ij}^{\mathrm{T}}{\dot{\mathit{q}}}_{n_{k}j}+{\dot{\mathit{q}}}_{ij}^{\mathrm{T}}{\ddot{\mathit{q}}}_{n_{k}j})+\\ \sum k_{n_k}({\dot{\mathit{q}}}_{ij}^T{\mathit{f}}_{a{n}_k})({\ddot{\mathit{q}}}_{ij}^{\mathrm{T}}{\mathit{f}}_{a{n}_k}+{\dot{\mathit{q}}}_{ij}^{\mathrm{T}}{\dot{\mathit{f}}}_{a{n}_k})=\\ -{\dot{\mathit{q}}}_{ij}^{\mathrm{T}}{C}_{ij}{\dot{\mathit{q}}}_{ij}-{\dot{\mathit{q}}}_{ij}^{\mathrm{T}}{K}_{ij}{\mathit{q}}_{ij}-{\dot{\mathit{q}}}_{ij}^{\mathrm{T}}{B}_{iuj}{\mathit{u}}_{ij}+\\ {\mathit{q}}_{ij}^{\mathrm{T}}{K}_{ij}{\dot{\mathit{q}}}_{ij}+{\mathit{q}}_{ij}^{\mathrm{T}}R{\dot{\mathit{q}}}_{ij}+\sum C_{i{n}_{k}j}({\dot{\mathit{q}}}_{ij}^T{\dot{\mathit{q}}}_{n_{k}j})·\\ ({\ddot{\mathit{q}}}_{ij}^{\mathrm{T}}{\dot{\mathit{q}}}_{n_{k}j}+{\dot{\mathit{q}}}_{ij}^{\mathrm{T}}{\ddot{\mathit{q}}}_{n_{k}j})+\sum k_{n_k}({\dot{\mathit{q}}}_{ij}^T{\mathit{f}}_{a{n}_k})·\\ ({\ddot{\mathit{q}}}_{ij}^{\mathrm{T}}{\mathit{f}}_{a{n}_k}+{\dot{\mathit{q}}}_{ij}^{\mathrm{T}}{\dot{\mathit{f}}}_{a{n}_k})=-{\dot{\mathit{q}}}_{ij}^{\mathrm{T}}{C}_{ij}{\dot{\mathit{q}}}_{ij}+\\ {\dot{\mathit{q}}}_{ij}^{\mathrm{T}}\alpha {\dot{\mathit{q}}}_{ij}+{\dot{\mathit{q}}}_{ij}^{\mathrm{T}}\beta {\dot{\mathit{q}}}_{ij}+{\mathit{q}}_{ij}^{\mathrm{T}}R{\dot{\mathit{q}}}_{ij}=\\ {\dot{\mathit{q}}}_{ij}^{\mathrm{T}}(\alpha +R){\mathit{q}}_{ij}+{\dot{\mathit{q}}}_{ij}^{\mathrm{T}}(\beta -C_{ij}){\dot{\mathit{q}}}_{ij}\end{array} （17）

由控制律的设置条件 可知Lyapunov函数对时间的一阶导数：

{\mathit{V}}^{\mathrm{\text{'}}} = {\dot{\mathit{q}}}_{ij}^{\mathrm{T}} ( \alpha + R ) {\mathit{q}}_{ij} + {\dot{\mathit{q}}}_{ij}^{\mathrm{T}} ( \beta - C_{ij} ) {\dot{\mathit{q}}}_{ij} （18）

恒小于零，根据Lyapunov稳定性理论，所设计的控制律能够保证系统的稳定性。

4 仿真分析

采用NX12.0建立卫星单侧大挠性太阳能帆板的虚拟样机模型，对其进行分析。设定帆板的材料各处均匀，忽略帆板与卫星中心刚体相连接的桁架部件，认为帆板在轨展开后是刚性锁定，帆板整体可以简化为板模型，作动器和传感器的安装位置相同，且控制器及其电路已经集成在帆板结构中，划分的子模块几何、材料参数相同。根据某项目需求及上文所求出的作动器、传感器最优安装位置，取主要的3阶模态进行振动抑制的仿真分析实验，模态振型如图7所示。具体的帆板固有频率如下所示：第1阶、第2阶、第3阶模态的固有频率分别为5.149 77、10.2914、14.2687 Hz。根据控制系统的分析，在Simulink中建立分布式控制方法的仿真模型并进行仿真验证。

4.1 外界持续干扰情况

为模拟卫星在轨工作的复杂环境，设置一个持续存在的外界干扰，该干扰为服从正态分布的随机信号，方差为0.05，采样时间为0.01 s，在帆板8个子模块处的干扰各不相同。由于子模块1和子模块5距离卫星的中心刚体最近，由其引起的振动对卫星的影响最大，所以仿真分析以子模块1和5为主，仿真时间为10 s。

首先，建立如式（12）所示的分布式控制器Simulink模型，所设计的分布式控制方法作用于大挠性太阳能帆板上，对帆板的模型进行仿真。由于帆板在5 s后的时间内保持稳定状态，故仅展示前5 s内的振动抑制控制图像。得到图8所示的子模块振动模态坐标收敛图像和图9所示的子模块输出控制力随时间变化图像。由于外界干扰的持续存在，所以最终振动抑制效果不可能为零，因此设定振动抑制的指标为振幅抑制达到97%。

为了验证本文分布式控制方法的有效性和优越性，将其与基于Lyapunov稳定性原理设计的集中式控制方法^{［19，20］}和基于LQR控制器的分布式控制方法^［10］进行对比分析，得到子模块1和子模块5的振动抑制仿真对比结果，如图10所示。以第1阶主模态为例，分别对比子模块1和子模块5的振动模态坐标收敛图像。从图10（a）中可以看出：对于太阳能帆板的第1阶模态，与传统的集中式控制方法和基于LQR的分布式控制方法相比，本文控制方法在振动抑制过程中，模态坐标的收敛时间较短，而且收敛过程中的模态坐标振动幅值较小。如表3所示，在持续外界干扰的情况下，振动幅值抑制达到97%的时间分别为：子模块1的时间为0.8928 s；子模块5的时间为0.8060 s；平均振动抑制时间为0.7812 s。本文方法抑振时间相比于基于LQR的分布式控制方法缩短了40%，相比于集中式控制方法振动抑制时间缩短了50%。由图10（b）（c）可以看出：对于太阳能帆板的第2、3阶模态，本文控制方法与其他两种控制方法相比，模态坐标收敛时间缩短，且抑振过程中的震荡减小，具有更好的振动抑制效果。

4.2 单次强冲击外界干扰情况

卫星大挠性帆板振动抑制最主要的目标是使距离中心刚体最近的子模块保持稳定，以保证附件振动对于中心刚体的运动影响最小。为了应对卫星在轨工作环境中的突发外界干扰情况，在已有持续外界干扰的基础上，对1号子模块添加单次强冲击外界干扰。为了模拟帆板在轨运行所处的状态，模拟受到短时大幅干扰，设定该干扰持续时间为0.1 s，沿z轴方向，加速度为5 m/s²：

式中： T_{\mathrm{d}}为短时大幅干扰。

将本文设计的分布式控制方法应用在太阳能帆板模型上，进行Simulink仿真，同样与基于LQR的分布式控制方法和集中式控制方法进行仿真对比分析，可以得到如图11所示的子模块振动模态坐标收敛图像和如图12所示的子模块输出控制力变化图像，以及如图13所示的以前3阶模态的模态坐标收敛图像。

由图11可知：由于5 s时在子模块1上施加了一个单次强冲击外界干扰，对于太阳能帆板的第1阶主模态影响较大，产生了一个较大的尖峰震荡，并且在其邻接的子模块5上第5 s时也出现了一个较小的尖峰振荡，这是由于邻接耦合影响的原因。由于该干扰是沿z轴方向，所以其对于第1阶主模态影响较大。

图12的子模块输出控制力同样由于设置干扰为沿z轴方向，对第2、3阶模态影响较小。振动抑制的细节如图13所示，图13（a）的第1阶模态振动在0.6 s内得到抑制，帆板的模态坐标振动幅值抑制达到97%。由图13（b）（c）的第2、3阶模态振动抑制结果可以看出：本文控制方法与其他两种控制方法相比，模态坐标收敛的时间缩短且抑振过程中的震荡减小，具有更好的振动抑制效果。

本文分布式控制方法在遇到单次强冲击外界干扰时具有更快的反应速度，产生的干扰尖峰较小。如表4所示，对于子模块1而言，振动抑制时间相比基于LQR的分布式控制方法减少了0.22 s，相比集中式控制方法减少了0.37 s。

5 结束语

本文针对卫星大挠性太阳能帆板的分布式振动控制方法进行了研究。按作动器、传感器在帆板上的安装位置进行子模块划分，基于有限元离散化分析，并考虑邻接子模块振动产生的耦合影响，建立了子模块的动力学模型。在此基础上，基于Lyapunov稳定性理论，设计了控制器，并添加协同项形成分布式振动抑制控制方法。仿真结果表明：①本文分布式协同控制方法能够有效实现卫星大挠性帆板在轨工作时的振动抑制。②本文控制方法相比于传统基于LQR的分布式控制方法和集中式控制方法，在同样的外界干扰下具有更短的稳定收敛时间。在持续外界干扰情况下，抑振时间分别缩短了40%和50%；在单次强冲击干扰情况下，抑振时间分别缩短了27%和37%。具有更良好的稳定控制效果。③分布式控制方法具有良好的鲁棒性，在单次强冲击的外界干扰下，也能维持系统的稳定，抑振时间缩短了25%。

参考文献

[1] Nadafi R, Kabganian M, Kamali A, et al. Super-twisting sliding mode control design based on Lyapunov criteria for attitude tracking control and vibration suppression of a flexible spacecraft[J]. Measurement and Control, 2019, 52(7/8): 814-831.

[2] Liu Feng, Yue Bao-zeng, Zhao Liang-yu. Attitude dynamics and control of spacecraft with a partially filled liquid tank and flexible panels[J]. Acta Astronautica, 2018, 143: 327-336.

[3] Yuan Q, Liu Y, Qi N. Active vibration suppression for maneuvering spacecraft with high flexible appendages[J]. Acta Astronautica, 2017, 139: 512-520.

[4] Luo Y J, Xu M L, Yan B, et al. PD control for vibration attenuation in Hoop truss structure based on a novel piezoelectric bending actuator[J]. Journal of Sound & Vibration, 2015, 339: 11-24.

[5] Wang Z, Xu M, Jia Y, et al. Vibration suppression-based attitude control for flexible spacecraft[J]. Aerospace Science & Technology, 2017, 70: 487-496.

[6] Rahman N U, Alam M N, Ansari J A. An experimental study on dynamic analysis and active vibration control of smart laminated plates[J]. Materials Today: Proceedings, 2021, 46: 9550-9554.

[7] Tian J, Guo Q, Shi G. Laminated piezoelectric beam element for dynamic analysis of piezolaminated smart beams and GA-based LQR active vibration control[J]. Composite Structures, 2020, 252: No.112480.

[8] Hu Q, Ma G, Li C. Active vibration control of a flexible plate structure using LMI-based H_∞ output feedback control law[C]//Fifth World Congress on Intelligent Control and Automation, Hangzhou, China, 2004: No.8369247.

[9] 苗双全, 丛炳龙, 刘向东. 基于输入成形的挠性航天器自适应滑模控制[J]. 航空学报, 2013, 34(8): 1906-1914.

[10] Miao Shuang-quan, Cong Bing-long, Liu Xiang-dong. Adaptive sliding mode control of flexible spacecraft on input shaping[J]. Acta Aeronautica Et Astronautica Sinica, 2013, 34(8): 1906-1914.

[11] Wang E, Wu S, Liu Y, et al. Distributed vibration control of a large solar power satellite[J]. Astrodynamics, 2019, 3(2): 189-203.

[12] Ji N, Liu J. Distributed vibration control for flexible spacecraft with distributed disturbance and actuator fault[J]. Journal of Sound and Vibration, 2020, 475: No.5274.

[13] Li Q, Yang H, Zhao D, et al. Fault-tolerant control and vibration suppression of flexible spacecraft: An interconnected system approach[J]. Chinese Journal of Aeronautics, 2020, 33(7): 2014-2023.

[14] Nakka Y, Chung S J, Allison J, et al. Nonlinear attitude control of a spacecraft with distributed actuation of solar arrays[J]. Journal of Guidance, Control, and Dynamics, 2019, 42(3): 458-475.

[15] Chen T, Shan J, Wen H. Distributed passivity-based control for multiple flexible spacecraft with attitude-only measurements[J]. Aerospace Science and Technology, 2019, 94: No.105408.

[16] 韩泽强. 大挠性卫星高精度控制关键技术研究[D]. 哈尔滨:哈尔滨工业大学航天学院,2020.

[17] Han Ze-qiang. Research on key technologies of high-precision satellite high-precision control[D]. Harbin: School of Astronautics, Harbin Institute of Technology, 2020.

[18] Muhammad A K, Wang X G, Cui N G, et al. A criterion for optimal sensor placement for minimizing spillover effects on optimal controllers[J]. Journal of Vibration and Control, 2018, 24(8): 1469-1487.

[19] 刘潇翔, 石恒, 王思野. 挠性空间结构的密集模态特性及影响分析[J].空间控制技术与应用,2017,43(01):11-16.

[20] Liu Xiao-xiang, Shi Heng, Wang Si-ye. An analysis on characteristics and impacts of close modes in flexible space structures[J]. Aerospace Control and Application, 2017, 43(1): 11-16.

[21] Åström K J, Murray R M. Feedback Systems: An Introduction for Scientists and Engineers[M]. Princeton: Princeton University Press, 2010.

[22] He W, Ge S Z S. Dynamic modeling and vibration control of a flexible satellite[J]. IEEE Transactions on Aerospace & Electronic Systems, 2015, 51(2): 1422-1431.

[23] 陆栋宁, 刘一武. 带旋转挠性太阳帆板卫星自适应控制[J]. 航天控制, 2014, 32(1): 49-54.

[24] Lu Dong-ning, Liu Yi-wu. Adaptive control of the spacecraft with a rotating flexible solar array[J]. Aerospace Control, 2014, 32(1): 49-54.

基金资助

国家自然科学基金项目(62073150)

国家自然科学基金项目(62073151)

国家自然科学基金项目(62001185)

吉林省科技厅项目(20210101177JC)