-
摘要: 面向空间攻防等任务的航天器通常安装微波、激光等大功率对抗载荷, 未来航天器需要装备大型挠性太阳能帆板. 针对挠性航天器姿态机动过程中存在外部干扰、执行机构饱和及挠性附件振动且挠性模态不易直接测量等问题, 提出带挠性附件航天器的全驱姿态控制方法. 首先, 建立挠性航天器全驱姿态控制模型. 其次, 基于扩展非线性观测器(Extended nonlinearity observer, ENO)与努斯鲍姆增益调节设计一种抗饱和的姿态控制鲁棒算法. 将外部扰动、挠性振动和输入饱和函数饱和估计误差作为复合干扰, 采用非线性干扰观测器对其进行有效补偿. 在直接参数设计线性控制参数基础上, 扩展非线性观测器负责对挠性航天器产生的挠性振动等非线性进行实时估计和补偿, 努斯鲍姆函数辅助控制器输出力矩避免饱和, 并利用李雅普诺夫方法严格证明闭环系统的稳定性. 最后通过数学仿真验证该方法不仅能够实现执行机构饱和约束条件下的姿态控制, 还能有效抑制挠性结构的振动, 为探索未来带有大型挠性附件航天器姿态控制新的方法提供参考.Abstract: For space missions such as space attack and defense, spacecrafts are usually equipped with high-power counter loads such as microwaves and lasers, thus future spacecrafts will need to be equipped with large-scale flexible solar panels. To solve the problems of external disturbances, actuator saturation, and vibrations of flexible attachments difficult to measure directly during attitude maneuver of flexible spacecraft, a fully actuated-based attitude control method for flexible spacecraft is proposed. First, a fully actuated attitude model for flexible spacecraft is established. Second, a robust attitude control algorithm against saturation is designed based on the extended nonlinear observer (ENO) and Nussbaum gain adjustment. The external disturbances, flexible vibrations, and the approximation error of the input saturation function are considered composite disturbances, which are efficiently compensated by the nonlinear disturbance observer. Based on the direct parametric design of linear control parameters, the ENO is responsible for the real-time estimation and compensation of nonlinearities such as flexible vibrations generated by the flexible spacecraft. And the Nussbaum function is used for the adjustment of the controller output torque magnitude to avoid actuator saturation. The stability of the closed-loop system is rigorously demonstrated using the Lyapunov method. Finally, it is verified by mathematical simulation that this method can not only realize the attitude control under the actuator saturation constraint but also effectively suppress the flexible vibrations, which contributes to exploring a novel method for the attitude control of spacecraft with large-scale flexural attachments in the future.
-
近年来, 针对 “星链”、“一网”等为典型代表的低轨大规模星座的监测、对抗[1]已成为未来争夺天基信息主导权的迫切需求. 考虑到微波、激光等对抗载荷的大功率要求, 未来航天器通常需要安装大型太阳能帆板. 大挠性、低阻尼附件与中心刚体之间存在相互耦合作用, 其产生的挠性振动会降低航天器姿态指向精度; 另一方面, 航天器在轨运行过程中不可避免地会受到各种复杂环境干扰的影响, 而且有效载荷运动、燃料消耗等因素会导致航天器转动惯量的不确定性. 此外, 受执行器的体积及功耗约束, 其存在输出饱和问题. 如果在控制设计中忽略饱和影响, 往往会降低控制系统的性能或导致不稳定. 因此, 设计一种鲁棒饱和姿态控制方法来完成挠性航天器的姿态机动控制和振动抑制具有重要意义.
在过去十几年中, 挠性航天器的姿态控制涌现出大量方法, 包括反馈控制[2-3]、自适应控制[4]、变结构控制[4-5]、容错控制[6-7]和扰动补偿控制[8-9]. 在实际应用中, 几乎所有控制执行机构都受到幅度和速率限制, 忽略这种限制的控制设计往往会导致不符合实际系统的瞬态响应, 降低闭环性能, 甚至可能导致闭环不稳定. 因此, 围绕饱和约束的控制问题是一个重要课题, 并已经有许多解决方法[8-10].
针对非线性系统, 除上述传统方法外, 计算力矩法也是一种功能强大的非线性控制器[11-12], 其广泛应用于机器人机械臂的控制. 该方法以反馈线性化为基础, 当模型动态参数精确已知时, 该方法控制性能优良, 但当模型存在参数不确定性、外部干扰时, 该控制器的性能锐减.
以上方法大都属于状态空间方法, 挠性航天器姿态动力学模型是一个经典的二阶欠驱系统, 传统的状态空间方法是将二阶系统转化为一阶系统, 再用现有的一阶方法进行处理. 但刚体姿态和柔性模态对应的高维状态空间模型计算困难, 欠驱动系统的非线性特性无法充分捕捉, 并且控制律设计繁琐. 全驱系统(Fully actuated system, FAS)理论正是在这样的背景下产生的. 针对状态空间法的问题, 段广仁院士于2020年首次提出FAS[13]. 与状态空间法不同, FAS理论保证了在物理环境不变的情况下, 无论非线性项的复杂性如何, 均能构造出数学意义上的全驱动结构来补偿系统的动态特性. 另外, 不同于计算力矩法对动力学模型高精度要求, FAS方法可对不确定模型、存在外部干扰的系统进行控制, 实现对系统不确定性、外部干扰等强鲁棒性控制.
自FAS诞生以来, 众多学者对FAS控制理论和应用进行了大量研究. 段广仁院士在FAS方法框架下研究了鲁棒控制[14]、自适应控制[15]、鲁棒自适应控制[16]等. 随着理论的发展和深入, FAS开始在工程中得到应用. 段广仁院士开发了一种航天器姿态稳定最优控制器[17]. 此外, FAS还应用于充液航天器的姿态控制[18]、刚性航天器[19]的容错跟踪控制. 只要从理论上选择适当的控制参数, 基于FAS的控制器就可以获得任意质量的动态性能. 但由于设计自由度大, 导致设计参数选取较为困难.
构造挠性航天器数学意义上的全驱系统模型不仅对深入理解挠性航天器的非线性动态特性和控制机制有着启发意义, 而且对工程上其他诸如柔性臂机器人等欠驱系统控制的解决具有延伸意义, 为欠驱系统控制研究提供参考. 鉴于执行机构输出力矩限幅与动态特性对闭环系统的瞬态响应的影响, 探索在控制输入饱和约束下的基于FAS理论的挠性航天器姿态控制方法具有重要意义. 然而, 目前此方面研究仍相对较少. 因此, 有必要基于FAS理论对在控制输入饱和约束下的挠性航天器姿态控制方法进行研究.
鉴于以上分析, 本文针对挠性航天器姿态机动 和附件振动抑制问题, 综合考虑外界干扰、执行机构饱和及挠性附件振动及振动模态不易直接测量等问题, 从全驱系统理论出发, 构建考虑执行器饱和、外部干扰等的数学意义上的挠性航天器全驱姿态模型, 基于扩展非线性观测器(Extended nonlinearity observer, ENO)与努斯鲍姆增益调节技术, 提出一种抗饱和的姿态控制与振动抑制参数控制方法. 该方法基于FAS方法, 采用直接参数法辅助反馈控制, 进行系统设计自由度的挖掘与控制参数选择, 一定程度上缓解了上述控制参数难以调节问题; 设计扩展非线性状态观测器观测包括系统惯性耦合项、挠性振动、外部干扰等非线性, 并基于该观测器设计非线性补偿前馈控制部分; 在输入饱和约束下, 采用努斯鲍姆增益调节技术对控制力矩进行处理, 最大程度降低输入饱和影响. 数学仿真验证所提控制器的有效性, 即该控制器能有效地抑制外部干扰、挠性结构的振动等, 在实现姿态机动的同时, 保证其控制输出满足执行机构饱和约束要求, 具有较好的控制性能.
1. 数学模型及控制问题描述
1.1 姿态运动学与动力学模型
假设挠性航天器所需发电功率为1 kW, 根据太阳能电池板发电效率约140 ~ 150 $ \text{W}/\text{m}^{2} $, 其简易模型的几何参数设计如图1, 卫星本体及太阳能帆板均采用蜂窝铝板, 支撑架采用碳纤维材料.
本文采用修正罗德里格斯参数(Modified rodrigues parameters, MRPs) $ {\boldsymbol{\sigma}}{\in{\bf{R}}^3} $来描述航天器姿态, 航天器的运动学方程为:
$$ {\boldsymbol{\dot{\sigma}}}={{{\boldsymbol{G}}}}({{{\boldsymbol{\sigma}}}}){{{\boldsymbol{\omega}}}} $$ (1) 式中, $ {\boldsymbol{G}}({\boldsymbol{\sigma}})\;=\;0.25[(1\;-\;{\boldsymbol{\sigma}}^{\text{T}}{\boldsymbol{\sigma}}){\boldsymbol{I}}_{\text{3}}\;+\;2{\boldsymbol{\sigma}}^{\times}\;+\;2({\boldsymbol{\sigma}}\;\times {\boldsymbol{\sigma}})^{\text{T}}] \in {\bf{R}}^{3\times3} $为坐标转换矩阵, $ {\boldsymbol{\omega}}=[{\omega_{x}},\;{\omega_{y}},\;{\omega_{z}}]^{\text{T}}\in {\bf{R}}^{3} $ 为航天器相对于惯性系在本体系下的姿态角速度. $ {\boldsymbol{I}}_{{3}} $为单位阵, 斜对称阵$ \boldsymbol{\sigma}^{\times}=[0,-\sigma_z,\; \sigma_y;\; \sigma_z,\; 0, -\sigma_x;\; -\sigma_y,\; \sigma_x,\;0] $.
鉴于航天器诸如太阳能帆板、天线等挠性附件带来的结构振动, 其姿态动力学模型表示为:
$$ \begin{cases} {\boldsymbol{J}}\dot{{\boldsymbol{\omega}}}+{\boldsymbol{\delta}}^\text{T}\ddot{{\boldsymbol{\eta}}}+{\boldsymbol{\omega}}^\times({\boldsymbol{J}}{\boldsymbol{\omega}}+{\boldsymbol{\delta}}^\text{T}\dot{{\boldsymbol{\eta}}})={\boldsymbol{u}}+{\boldsymbol{T}}_d\\ \ddot{{\boldsymbol{\eta}}}+{\boldsymbol{C}}\dot{{\boldsymbol{\eta}}}+{\boldsymbol{K}}{\boldsymbol{\eta}}=-{\boldsymbol{\delta}}\dot{{\boldsymbol{\omega}}} \end{cases} $$ (2) 式中, $ {\boldsymbol{J}}\in{\bf{R}}^{3\times3} $为航天器转动惯量矩阵; $ {\boldsymbol{\delta}}\in{\bf{R}}^{4\times3} $为中心刚体与挠性附件的耦合矩阵; $ {\boldsymbol{\eta}}\in{\bf{R}}^{4\times1} $为挠性振动模态; $ {\boldsymbol{u}}\in{\bf{R}}^{3\times1} $为控制力矩; $ {\boldsymbol{T}}_d\in{\bf{R}}^{3\times1} $为外部干扰力矩; 阻尼矩阵$ {\boldsymbol{C}} $、刚度矩阵$ {\boldsymbol{K}} $分别为
$$ \begin{cases}{\boldsymbol{C}} = {\mathrm{diag}}\left\{ {2\xi {\omega _{n1}},{\rm{ }}2\xi {\omega _{n2}},{\rm{ }}2\xi {\omega _{n3}},{\rm{ }}2\xi {\omega _{n4}}} \right\} \in {{\bf{R}}^{{\rm{n}} \times {\rm{n}}}}\\{\boldsymbol{K}} = {\mathrm{diag}}\left\{ {\omega _{n1}^2,{\rm{ }}\omega _{n2}^2,{\rm{ }}\omega _{n3}^2,{\rm{ }}\omega _{n4}^2} \right\} \in {{\bf{R}}^{{\rm{n}} \times {\rm{n}}}}\end{cases} $$ (3) 其中, $ {\xi_i} $、$ \omega_{ni} $分别为航天器$ i\left(i=1,\;2,\;3,\;4\right) $阶阻尼比、模态固有频率, $ n $为模态阶数.
考虑飞轮等执行器饱和, 式(2)中真实的系统控制力矩$ {\boldsymbol{u}} $一般表示为$ {\boldsymbol{u}}({\boldsymbol{\nu}})=[\text{sat}(\nu_1),\;\text{sat}(\nu_2), \text{sat}(\nu_3)]^\text{T} $, $ {\boldsymbol{\nu}}\;\in\;{\bf{R}}^{3\times1} $ 为控制指令力矩, $ \text{sat}(\nu_i)\;= \text{sgn}(\nu_i)\cdot \min\{\left|\nu_i\right|,\;u_{i,\;\max}\} $, $ {u}_{i,\;\max} $为对应执行器最大输出力矩. 本文选用双曲正切函数$ {\boldsymbol{g}}({\boldsymbol{\nu}})\in{\bf{R}}^{3\times1} $近似饱和非线性, 因此, 实际控制力矩$ {\boldsymbol{u}}({\boldsymbol{\nu}}) $表示为:
$$ {{{\boldsymbol{u}}}}({{{\boldsymbol{\nu}}}})={{{\boldsymbol{g}}}}({{{\boldsymbol{\nu}}}})+{{{\boldsymbol{e}}}}_s({{{\boldsymbol{\nu}}}}) $$ (4) 其中, $ {g}_i(\nu_i)=u_{i,\;\max}\times\tanh(\nu_i/u_{i,\;\max}),\;i=1,\;2,\;3 $; 因为 $ \left|e(\nu_i)\right|=\left|\operatorname{sat}(\nu_i)-{ g}_i(\nu_i)\right|\leq u_{i,\;\max}(1- \tanh(1)), i=1,\;2,\;3 $, 近似误差向量$ {\boldsymbol{e}}_s({\boldsymbol{\nu}}) $有界. 另外, 定义对角阵$ {\boldsymbol{H}}\in{\bf{R}}^{3\times3} $, 其对角元$ h_i={g}_i(\nu_i)/\nu_i,\; h_i\in (0,\;1] $. 因此, 考虑执行器饱和, 实际控制力矩(4)可转换为
$$ {\boldsymbol{u}}({\boldsymbol{\nu}})={\boldsymbol{H}}{\boldsymbol{\nu}}+{\boldsymbol{e}}_s({\boldsymbol{\nu}}) $$ (5) 注 1. 为避免奇异, 当$ \nu_{i}=0 $时, 令$ h_i=1 $.
因此, 考虑输入饱和, 挠性航天器姿态动力学模型(2)可表示为:
$$ \begin{cases}{\boldsymbol{J}}\dot{{\boldsymbol{\omega}}}+{\boldsymbol{\delta}}^\text{T}\ddot{{\boldsymbol{\eta}}}+{\boldsymbol{\omega}}^\times({\boldsymbol{J}}{\boldsymbol{\omega}}+{\boldsymbol{\delta}}^\text{T}\dot{{\boldsymbol{\eta}}})={\boldsymbol{H}}{\boldsymbol{\nu}}+{\boldsymbol{d}}\\\ddot{{\boldsymbol{\eta}}}+{\boldsymbol{C}}\dot{{\boldsymbol{\eta}}}+{\boldsymbol{K}}{\boldsymbol{\eta}}=-{\boldsymbol{\delta}}\dot{{\boldsymbol{\omega}}}\end{cases} $$ (6) 式中, $ {\boldsymbol{d}}={\boldsymbol{e}}_s({\boldsymbol{\nu}})+{\boldsymbol{T}_d} $为包括外部干扰、饱和函数饱和估计误差的总扰动. 将航天器姿态模型(1)与(6)转化为二阶全驱系统(Second order fully actuated system, SOFAS)为
$$ \begin{cases} \dot{{\boldsymbol{\sigma}}}={\boldsymbol{f}}+{\boldsymbol{Bu}}\\{\boldsymbol{f}}={\boldsymbol{f}}_{1}+{\boldsymbol{f}}_{2}+{\boldsymbol{f}}_{3}\\ {\boldsymbol{f}}_{1}=[\dot{{\boldsymbol{G}}}{\boldsymbol{G}}^{-1}-{\boldsymbol{G}}({\boldsymbol{J}}-{\boldsymbol{\delta}}^{\text{T}}{\boldsymbol{\delta}})^{-1}{\boldsymbol{\omega}}^\times{\boldsymbol{J}}{\boldsymbol{G}}^{-1}]\dot{{\boldsymbol{\sigma}}}\\ {\boldsymbol{f}}_{2}={\boldsymbol{G}}({\boldsymbol{J}}-{\boldsymbol{\delta}}^{\text{T}}{\boldsymbol{\delta}})^{-1}({\boldsymbol{\delta}}^{\text{T}}{\boldsymbol{C}}\dot{{\boldsymbol{\eta}}}+{\boldsymbol{\delta}}^{\text{T}}{\boldsymbol{K}}{\boldsymbol{\eta}}-{\boldsymbol{\omega}}^\times{\boldsymbol{\delta}}^{\text{T}}\dot{{\boldsymbol{\eta}}})\\ {\boldsymbol{f}}_{3}={\boldsymbol{G}}({\boldsymbol{J}}-{\boldsymbol{\delta}}^{\text{T}}{\boldsymbol{\delta}})^{-1}{\boldsymbol{d}}={\boldsymbol{G}}({\boldsymbol{J}}-{\boldsymbol{\delta}}^{\text{T}}{\boldsymbol{\delta}})^{-1}({\boldsymbol{e}}_{s}+{\boldsymbol{T}}_{d})\\ {\boldsymbol{Bu}}={\boldsymbol{G}}({\boldsymbol{J}}-{\boldsymbol{\delta}}^{\text{T}}{\boldsymbol{\delta}})^{-1}{\boldsymbol{u}}={\boldsymbol{G}}({\boldsymbol{J}}-{\boldsymbol{\delta}}^{\text{T}}{\boldsymbol{\delta}})^{-1}{\boldsymbol{Hv}} \end{cases} $$ (7) 式中, $ {\boldsymbol{f}}_1、{\boldsymbol{f}}_2、{\boldsymbol{f}}_3\in{\bf{R}}^3 $分别表征刚体姿态运动惯性耦合、刚柔耦合引起的非线性振动、外部干扰与饱和估计误差. 另外, $ {\boldsymbol{f}}({\boldsymbol{\sigma}},\;\dot{{\boldsymbol{\sigma}}})\in{\bf{R}}^3 $、$ {\boldsymbol{B}}({\boldsymbol{\sigma}},\;\dot{{\boldsymbol{\sigma}}})\in{\bf{R}}^{3\times3} $为连续矢量函数, 由于雅可比矩阵$ {\boldsymbol{G}} $可逆, 转动惯量矩阵$ {\boldsymbol{J}} $为正定矩阵, 因此, $ {\boldsymbol{B}}\neq0 $, 即系统(7)满足全驱条件.
1.2 定义与引理
定义 1[20]. 任何连续偶函数$ N(\chi) $若满足下述条件, 称为努斯鲍姆型函数. 本文将其应用于考虑输入饱和的姿态控制问题可削弱执行机构饱和对控制器性能的不利影响.
$$ \left\{\begin{aligned} &\lim_{s\to\infty}\sup\frac1s\int_0^sN(\chi)\text{d}\chi=+\infty\\ &\lim_{s\to\infty}\inf\frac1s\int_0^sN(\chi)\text{d}\chi=-\infty \end{aligned}\right. $$ 典型的努斯鲍姆型函数有$ \chi^2\cos(\chi) $、$ \chi^2\sin(\chi) $、$ \text{e}^{\chi^2} $等. 本文选择的努斯鲍姆函数如式(8)所示[15-17].
$$ N(\chi)=\text{e}^{\frac{\chi^2}{2}}(\chi^2+2)\sin(\chi)+1 $$ (8) 引理 1[20-21]. 设$ V(t) $和$ \chi_i(t) $是定义在$ \left[0,\;t_f\right) $上的光滑函数, 且$ \chi_i(t)=0,\; N(\cdot) $是由式(8)定义的努斯鲍姆型函数. 若下述不等式成立,
$$ \begin{split} V(t)\leq\;& c_0+\text{e}^{-c_1}\sum_{i=1}^N\int_0^t(-g_i(\tau)N(\chi_i(\tau))\;+\\ &1)\dot{\chi}_i(\tau)\text{e}^{c_1\tau}\text{d}\tau \end{split} $$ 式中, $ c_{0} $为有界常数, $ c_{1}>0 $, 时变参数$ g_{i}(t) $满足$ g_{i}(t)=\left\{g_{i}(t)\Bigg|\begin{aligned}&g_{i}(t)\in\text{ 未知区间}I\\ &I=[g^{-},\;0)\cup(0,\;g^{+}] \end{aligned}\right\} $, 则$ V(t) $, $ \chi_{i}(t) $和$\sum\nolimits_{i=0}^N {\int_0^t {{g_i}(\tau )} } N({\chi_i}(\tau ))\cdot {\dot \chi _i}(\tau ){{\rm{d}}}\tau$在$ [0,\;t_f) $上有界.
引理 2[22]. 设$ \varPsi,\;V:[0,\;\infty)\to{\bf{R}} $, 且$ t_0\in(0,\;t) $, 若存在常数$ \alpha>0 $, 使得$ \dot{V}\leq-\alpha V+\Psi $, 则
$$ V(t)\leq\text{e}^{-\alpha(t-t_0)}V(t_0)+\int_0^t\text{e}^{-\alpha(t-\tau)}\Psi(\tau)\text{d}\tau $$ (9) 引理 3[23]. $ \forall {{\boldsymbol x},\;{\boldsymbol y}}\in{\bf{R}}^{n},\;{\boldsymbol{P}}>{\bf{0}} $, 矩阵$ {\boldsymbol{A}} $, $ {\boldsymbol{U}} $ 满足下述不等式
$$ 2{\boldsymbol{x}}^\text{T}{\boldsymbol{AUy}}\leq {\boldsymbol{x}}^\text{T}{\boldsymbol{APA}}^\text{T}{\boldsymbol{x}}+{\boldsymbol{y}}^\text{T}{\boldsymbol{U}}^\text{T}{\boldsymbol{P}}^{-1}{\boldsymbol{Uy}} $$ (10) 2. 全驱姿态控制器设计
为了方便后续控制器的设计, 作如下假设: 1)星体转动惯量矩阵为常矩阵; 2)外部环境扰动、模态位移和模态速度未知但有界且此类非线性干扰变化较为缓慢, 即变换后的总扰动$ {\boldsymbol{f}} $有界可微, 并且存在正标量$ \delta d_1 $和$ \delta d_2 $, 使得$ \left\|{\boldsymbol{f}}\right\|\leq\delta d_1 $与$ \|\dot{\boldsymbol f}\|\leq\delta d_2 $成立.
本文的控制目标可描述为: 针对存在外界干扰、控制受限、挠性附件振动及挠性模态不易直接测量情况下的挠性航天器姿态控制系统式(1)、(6), 在假设1 ~ 2条件下, 设计一类非线性饱和姿态控制器, 如图2所示.
控制器基于全驱系统理论, 结构上包括线性反馈主部与非线性补偿部分. 采用直接参数法设计线性反馈参数, 扩展非线性观测器估计未知非线性, 最终基于努斯鲍姆增益调节与饱和函数对控制输出进行抗饱和处理. 其在完成姿态机动的同时, 抑制挠性附件的结构振动, 并保证控制输出满足执行机构饱和受限要求, 且对外部干扰具有较好的鲁棒性.
2.1 扩展非线性扰动观测器设计
已知系统(7)状态$ {\boldsymbol{\sigma}} $、$ \dot{{\boldsymbol{\sigma}}} $, 将其未知非线性部分设为扩展状态$ {\boldsymbol{z}}_{{3}}={\boldsymbol{f}}(\dot{{\boldsymbol{\sigma}}}) $, 定义$ {\boldsymbol{z}}_{{1}}={\boldsymbol{\sigma}},\;{\boldsymbol{z}}_{{2}}=\dot{{\boldsymbol{\sigma}}} $, 且由假设2, $ \dot{{\boldsymbol{z}}}_3={\boldsymbol{k}}(t) $未知有界. 因此, 式(7)可转化为
$$ \begin{cases} \dot{{\boldsymbol{z}}}_1={\boldsymbol{z}}_{{2}}\\\dot{{\boldsymbol{z}}}_2={\boldsymbol{z}}_3+{\boldsymbol{Bu}}\\\dot{{\boldsymbol{z}}}_3={\boldsymbol{k}}(t) \end{cases} $$ (11) 基于原有系统, 设计一个与之平行的ENO系统:
$$ \begin{cases}\dot{\hat{{\boldsymbol{z}}}}_1={\boldsymbol{\hat{z}}}_2+\beta_1{\boldsymbol{\varphi}}_1({\boldsymbol{e}}_1)\\\dot{\hat{{\boldsymbol{z}}}}_2=\hat{{\boldsymbol{z}}}_3+{\boldsymbol{Bu}}+\beta_2{\boldsymbol{\varphi}}_2({\boldsymbol{e}}_1)\\\dot{\hat{{\boldsymbol{z}}}}_3=\beta_3{\boldsymbol{\varphi}}_3({\boldsymbol{e}}_1) \end{cases} $$ (12) 式中, 定义${{\boldsymbol{\varphi}} _i}({{\boldsymbol{e}}_1}) = {[{\varphi _i}({e_{1x}}),\;{\varphi _i}({e_{1y}}),\;{\varphi _i}({e_{1z}})]^{\mathrm{T}}}, \;i= 1,\;2,\;3,\;4$. $ {\boldsymbol{\hat{z}}}_1 $、$ {\boldsymbol{\hat{z}}}_2 $、$ {\boldsymbol{\hat{z}}}_3 $分别为$ {\boldsymbol{z}}_1 $、$ {\boldsymbol{z}}_2 $、$ {\boldsymbol{z}}_3 $的估计, $\varphi_{{i}}(\cdot)$为修正误差指数增益函数, 其定义如下
$$ \varphi_{i}(e,\;\alpha_i,\;\delta_i)=\begin{cases}\left|\text{e}\right|^{\alpha_i}\text{sgn}(e),&\left|e\right|>\delta_i\\\dfrac{e}{\delta_i^{1-\alpha_i}},&\left|e\right|\leq\delta_i\end{cases} $$ 其中, e为估计误差, $ {\alpha_i}\in(0,\;1) $, 保证在估计误差较小时$\varphi_{{i}}(\cdot)$具有高增益; $ \delta_i $是一个微小的正数, 用于限制原点邻域内的$\varphi_{{i}}(\cdot)$增益.
由平行系统(11) ~ (12), 定义$ {\boldsymbol{e}}_1=\hat{{\boldsymbol{z}}}_1-{\boldsymbol{z}}_1 $, 可得观测器误差动力学方程
$$ \begin{cases}\dot{{\boldsymbol{e}}}_1={\boldsymbol{e}}_2+\beta_1{\boldsymbol{\varphi}}_1({\boldsymbol{e}}_1)\\\dot{{\boldsymbol{e}}}_2={\boldsymbol{e}}_3+{\boldsymbol{B}}{\boldsymbol{u}}+\beta_2{\boldsymbol{\varphi}}_2({\boldsymbol{e}}_1)\\\dot{{\boldsymbol{e}}}_3=\beta_3{\boldsymbol{\varphi}}_3({\boldsymbol{e}}_1)-{\boldsymbol{k}}(t)\end{cases} $$ (13) 鉴于二阶全驱系统(7), 确定合适参数$ \beta_1 $, $ \beta_2 $, $ \beta_3 $, 使得系统在存在饱和的情况下满足
$$ \lim_{x\to\infty}{\boldsymbol{e}}_3(t)=\lim_{x\to\infty}[\hat{{\boldsymbol{z}}}_3-{\boldsymbol{z}}_3]=0 $$ (14) 2.2 控制律设计
2.2.1 线性反馈直接参数设计法
基于全驱系统理论框架设计关于系统(7)的姿态控制指令如下:
$$ \left\{\begin{aligned}&{\boldsymbol{v}}={\boldsymbol{N}}_D(\chi){\boldsymbol{\tau}}_a\\&{\boldsymbol{\tau}}_a={\boldsymbol{u}}_f+{\boldsymbol{u}}_d\end{aligned}\right. $$ (15) 式中, $ {\boldsymbol{N}}_D(\chi)={{\mathrm{diag}}}\{N(\chi_1),\;N(\chi_2),\;N(\chi_3)\} $; $ {\boldsymbol{\tau}}_a\in {\bf{R}}^{3\times1} $为控制律输出力矩; $ {\boldsymbol{u}}_f $为线性状态反馈部分, 以获得期望线性闭环系统; $ {\boldsymbol{u}}_d $用于补偿系统非线性$ {\boldsymbol{f}} $.
$$ {\boldsymbol{u}}_f=-{\boldsymbol{B}}^{-1}\left({\boldsymbol{A}}_0{\boldsymbol{\sigma}}+{\boldsymbol{A}}_1\dot{{\boldsymbol{\sigma}}}\right)+{\boldsymbol{v}}_{ex} $$ (16) 其中, $ {\boldsymbol{A}}_0 $、$ {\boldsymbol{A}}_1 $为待设计反馈增益矩阵, $ {\boldsymbol{v}}_{ex} $为外部输入前馈控制信号.
$$ {\boldsymbol{u}}_d=-{\boldsymbol{B}}^{-1}\hat{{\boldsymbol{z}}}_3=-({\boldsymbol{G}}({\boldsymbol{J}}-{\boldsymbol{\delta}}^{\text{T}}{\boldsymbol{\delta}})^{-1})^{-1}\hat{{\boldsymbol{z}}}_3 $$ (17) 式中, 针对系统非线性项$ {\boldsymbol{f}} $具有模糊性、不确定性等特点, 在控制律中引入对非线性的估计状态$ \hat{{\boldsymbol{z}}}_3 $.
控制器(15) ~ (17)应用于SOFAS模型(7), 可得常线性闭环系统
$$ \ddot{{\boldsymbol{\sigma}}}+{\boldsymbol{A}}_1\dot{{\boldsymbol{\sigma}}}+{\boldsymbol{A}}_0{\boldsymbol{\sigma}}=-{\boldsymbol{e}}_3+{\boldsymbol{B}}{\boldsymbol{\nu}}_{ex} $$ 鉴于本文方法设有针对外部干扰等非线性补偿前馈控制部分$ {\boldsymbol{u}}_d $, 因此不予考虑外部输入$ {\boldsymbol{\nu}}_{ex} $, 当非线性估计误差$ {\boldsymbol{e}}_3\to{\bf{0}} $, 上式转化为状态空间形式
$$ \dot{{\boldsymbol\sigma}}^{(0\sim1)}={\text{Φ}}({\boldsymbol{A}}_{0\sim1}){\boldsymbol\sigma}^{(0\sim1)} $$ 式中, $ {\text{Φ}}({\boldsymbol{A}}_{0\sim1})=\begin{bmatrix}{\boldsymbol{O}}_{3\times3}&{\boldsymbol{I}}_{3\times3}\\-{\boldsymbol{A}}_0&-{\boldsymbol{A}}_1\end{bmatrix} $. ${\boldsymbol{O}}_{3\times3} $表示${3\times3} $的零矩阵, ${\boldsymbol{I}}_{3\times3} $表示${3\times3} $的单位矩阵. 此时问题转化成寻找增益矩阵$ {\boldsymbol{A}}_{0\sim1}\in{\bf{R}}^{3\times3} $, 使得$ {\text{Φ}}({\boldsymbol{A}}_{0\sim1})\in {\bf{R}}^{6\times6} $为稳定阵. 其中, 矩阵$ {\boldsymbol{A}}_{0} $、$ {\boldsymbol{A}}_{1} $分别为比例微分增益项, 用于调整系统的动态特性与稳态性能. 一般地, 随着比例增益矩阵$ {\boldsymbol{A}}_{0} $增大, 系统响应速度变快, 但易引起系统振荡, 系统超调量增加, 而微分增益项矩阵$ {\boldsymbol{A}}_{1} $相当于系统的阻尼项, 能有效抑制系统振荡, 减小系统超调. 在控制过程中需要均衡调节矩阵$ {\boldsymbol{A}}_{0} $、$ {\boldsymbol{A}}_{1} $保证系统动态性能. 选定Hurwitz矩阵$ {\boldsymbol{F}}\in{\bf{R}}^{6\times6} $ [24], 任意选定$ {\boldsymbol{Z}}\in{\bf{R}}^{3\times6} $满足$ \det{\boldsymbol{V}}({\boldsymbol{Z}}, {\boldsymbol{F}})\neq 0 $, $ {\boldsymbol{A}}_{0\sim1}={\boldsymbol{Z}}{\boldsymbol{F}}^2{\boldsymbol{V}}^{-1} $.
设计线性状态反馈控制器
$$ \begin{cases}{\boldsymbol{u}}_f=[{\boldsymbol{K}}_p,\; {\boldsymbol{K}}_d][{\boldsymbol{\sigma}},\; \dot{{\boldsymbol{\sigma}}}]^\text{T}={\boldsymbol{B}}^{-1}{\boldsymbol{Z}}{\boldsymbol{F}}^2{\boldsymbol{V}}^{-1}[{\boldsymbol{\sigma}},\; \dot{{\boldsymbol{\sigma}}}]^\text{T}\\{\boldsymbol{V}}=({\boldsymbol{Z}},\; {\boldsymbol{Z}}{\boldsymbol{F}})^\text{T}\end{cases} $$ (18) 在参数矩阵$ {\boldsymbol{F}}\in{\bf{R}}^{6\times6}\text{、}{\boldsymbol{Z}}\in{\bf{R}}^{3\times6} $作用下, 模型(7)转换成
$$ {\boldsymbol{\dot{X}}}={\boldsymbol{VFV}}^{-1}{\boldsymbol{X}} $$ (19) 其中, 增广变量$ {\boldsymbol{X=}}\; [{\boldsymbol{\sigma}}^\text{T},\;\dot{{\boldsymbol{\sigma}}}^\text{T}]^\text{T} $.
2.2.2 基于努斯鲍姆增益调节的输入饱和控制法
这里努斯鲍姆函数自适应变量$ {\boldsymbol{\chi}} $与系统状态、控制指令相关, 通过$ {\boldsymbol{\chi}} $的调节可有效协调系统性能与系统能耗, 其满足更新律:
$$ \begin{cases}{\boldsymbol{\chi}}=-{{\zeta }}{\boldsymbol{s}}_{D}{\boldsymbol{\tau}}_a,\;{\boldsymbol{\chi}}(0)={\bf{0}}\\{\boldsymbol{s}}_{D}={{\rm{diag}}} \{\mathit{s}_1,\;\mathit{s}_2,\;\mathit{s}_3\}\\{\boldsymbol{s}}={\boldsymbol{X}}^{\text{T}}{\boldsymbol{P}}\left[{\bf{0}},\;{\boldsymbol{B}}\right]^{\text{T}}\in{\bf{R}}^{1\times3}&\end{cases} $$ (20) 式中, $ \zeta $为$ {\boldsymbol{\chi}} $更新因子, $ {\boldsymbol{\tau}}_a $为全驱系统理论框架下的控制律输出, $ {\boldsymbol{P}}\in{\bf{R}}^{6\times6} $的正定性在后续讨论.
3. 稳定性分析
定理 1. 考虑姿态控制系统(2), 当存在标量$ \varepsilon>0 $, 正定矩阵$ {\boldsymbol{P}} $满足式(21)时, 在实时更新(20)的控制器(16) ~ (18)和扩展状态观测器(12)作用下, 确保MRP姿态参量及角速度$ {\boldsymbol{\sigma}}、\dot{{\boldsymbol{\sigma}}} $收敛到平衡点附近的一个小区域.
$$ {\boldsymbol{PA}}_c+{\boldsymbol{A}}_c^\text{T}{\boldsymbol{P}}+\varepsilon {\boldsymbol{P}}<0 $$ (21) 证明. 构建系统(7) Lyapunov函数
$$ V=\frac12{\boldsymbol{X}}^\text{T}{\boldsymbol{PX}} $$ (22) 其中, $ {\boldsymbol P} > {\boldsymbol 0} $是矩阵不等式(Linear matrix inequalities, LMI)的解.
对系统(7)的$ {V} $关于时间求一阶导
$$ \begin{split} \dot{V}=\;&{\boldsymbol{X}}^\text{T}{\boldsymbol{P}}\begin{bmatrix}\dot{{\boldsymbol{\sigma}}}\\{\boldsymbol{f}}+{\boldsymbol{B}}{\boldsymbol{u}}\end{bmatrix}=\\ &{\boldsymbol{X}}^\text{T}{\boldsymbol{P}}\begin{bmatrix}\dot{{\boldsymbol{\sigma}}}\\(\dot{{\boldsymbol{G}}}{\boldsymbol{G}}^{-1}{\boldsymbol-}{\boldsymbol{G}}{\boldsymbol{J}}^{-1}{\boldsymbol{\omega}}^\times{\boldsymbol{J}}{\boldsymbol{G}}^{-1})\dot{{\boldsymbol{\sigma}}}\end{bmatrix}+ \\ &{\boldsymbol{X}}^\text{T}{\boldsymbol{P}}{\begin{bmatrix}{\boldsymbol{0}}\\-{\boldsymbol{GJ}}^{-1}({\boldsymbol{\delta}}^\text{T}\ddot{{\boldsymbol{\eta}}}+{\boldsymbol{\omega}}^\times{\boldsymbol{\delta}}^\text{T}\dot{{\boldsymbol{\eta}}})+{\boldsymbol{GJ}}^{-1}{\boldsymbol{d}}\end{bmatrix}}+\\ &{\boldsymbol{X}}^\text{T}{\boldsymbol{P}}{\begin{bmatrix}{\boldsymbol{0}}\\{\boldsymbol{GJ}}^{-1}{\boldsymbol{H}}{\boldsymbol{\nu}}\end{bmatrix}}= \\ &{\boldsymbol{X}}^\text{T}{\boldsymbol{P}}{\begin{bmatrix}\dot{{\boldsymbol{\sigma}}}\\(\dot{{\boldsymbol{G}}}{\boldsymbol{G}}^{-1}{\boldsymbol{-G}}{\boldsymbol{J}}^{-1}{\boldsymbol{\omega}}^\times{\boldsymbol{J}}{\boldsymbol{G}}^{-1})\dot{{\boldsymbol{\sigma}}}\end{bmatrix}}+ \\ &{\boldsymbol{X}}^\text{T}{\boldsymbol{P}}{\begin{bmatrix}{\boldsymbol{0}}\\-{\boldsymbol{GJ}}^{-1}({\boldsymbol{\delta}}^\text{T}\ddot{{\boldsymbol{\eta}}}+{\boldsymbol{\omega}}^\times{\boldsymbol{\delta}}^\text{T}\dot{{\boldsymbol{\eta}}})+{\boldsymbol{GJ}}^{-1}{\boldsymbol{d}}\end{bmatrix}}+ \\ &{\boldsymbol{X}}^\text{T}{\boldsymbol{P}}{\begin{bmatrix}{\boldsymbol{0}}\\{\boldsymbol{G}}{\boldsymbol{J}}^{-1}{\boldsymbol{\tau}}_a\end{bmatrix}}+{\boldsymbol{X}}^\text{T}{\boldsymbol{P}}{\begin{bmatrix}{\boldsymbol{0}}\\{\boldsymbol{G}}{\boldsymbol{J}}^{-1}\end{bmatrix}}({\boldsymbol{H}}{\boldsymbol{\nu}}-{\boldsymbol{\tau}}_a) \end{split} $$ (23) 将式(16) ~ (18)中的基于直接参数法的控制律$ {\boldsymbol{\tau}}_{a} $代入式(23), $ \dot{V} $满足
$$ \begin{split}\dot{V}=\;&{\boldsymbol{X}}^\text{T}{\boldsymbol{P}}[{\boldsymbol{A}}_c{\boldsymbol{X}}+{\boldsymbol{f}}-\hat{{\boldsymbol{z}}}_3]+{\boldsymbol{s}}({\boldsymbol{H}}{\boldsymbol{N}}_D{\boldsymbol{\tau}}_a-{\boldsymbol{\tau}}_a)=\\&\dot{V}_1+\dot{V}_2\\[-1pt]\end{split} $$ (24) 其中, $ {\boldsymbol{A}}_c={\boldsymbol{VFV}}^{-1} $, $ {\boldsymbol{s}} $定义见式(20). 另外令$ {\boldsymbol{\gamma}}= {\boldsymbol{f}}- \hat{{\boldsymbol{z}}}_3,\;\ {\boldsymbol{E}}_c=[{\bf{0}},\;{\boldsymbol{I}}_3]^\text{T} $, 可得$ \dot{V}_{1} $
$$ \begin{split} \dot{V}_{1}=\;& {\boldsymbol{X}}^\text{T}{\boldsymbol{P[}}{\boldsymbol A}_c{\boldsymbol{X}}+{\boldsymbol{f}}-\hat{{\boldsymbol{z}}}_3]=\\ &\frac12{\boldsymbol{X}}^\text{T}[{\boldsymbol{PA}}_c+{\boldsymbol{A}}_c^\text{T}{\boldsymbol{P]X}}+{\boldsymbol{X}}^\text{T}{\boldsymbol{PE}}_c{\boldsymbol{\gamma}}\leq\\ &\frac12{\boldsymbol{X}}^\text{T}[{\boldsymbol{PA}}_c+{\boldsymbol{A}}_c^\text{T}{\boldsymbol{P}}]{\boldsymbol{X}}+\frac\varepsilon2{\boldsymbol{X}}^\text{T}{\boldsymbol{PX}}\;+\\ &\frac1{2\varepsilon}{\boldsymbol{\gamma}}^\text{T}{\boldsymbol{E}}_c^\text{T}{\boldsymbol{PE}}_c{\boldsymbol{\gamma }} \end{split} $$ (25) 式中, 最后一个不等式由引理3推知, 其中$ {\boldsymbol{U}}= {\boldsymbol{E}}_{c} $, $ {\boldsymbol{x}}={\boldsymbol{X}},\;\ {\boldsymbol{y}}={\boldsymbol{\gamma }} $.
另外, 考虑到$ {\boldsymbol{\chi}} $的更新规律(20), 可得$ \dot{V}_{2} $
$$ \begin{split} \dot{V}_{2}=\;& {\boldsymbol{sHN}}_D(\chi){\boldsymbol{\tau}}_a-{\boldsymbol{\tau}}_a=\\ &{\boldsymbol{s}}({\boldsymbol{HN}}_D(\chi)-{\boldsymbol{I}}_3){\boldsymbol{\tau}}_a=\\ &-\frac1\zeta{\boldsymbol{s}}({\boldsymbol{HN}}_D(\chi)-{\boldsymbol{I}}_3){\boldsymbol{s}}_D^{-1}\dot{{\boldsymbol{\chi}}}=\\ &-\frac1\zeta\sum_{i=1}^3{\boldsymbol{s}}(h_i{\boldsymbol{N}}(\chi_i)-1){\boldsymbol{s}}_D^{-1}\dot{\chi}_i \end{split} $$ (26) 由式(25)与式(26), 式(24)可转化成
$$ \begin{split} \dot{V}=\;& \frac12{\boldsymbol{X}}^\text{T}[{\boldsymbol{PA}}_c+{\boldsymbol{A}}_c^\text{T}{\boldsymbol{P}}+\varepsilon {\boldsymbol{P}}]{\boldsymbol{X}}\;+ \\ &\frac12{\boldsymbol{X}}^\text{T}[{\boldsymbol{PA}}_c+{\boldsymbol{A}}_c^\text{T}{\boldsymbol{P}}]{\boldsymbol{X}}+\frac1{2{{\varepsilon}}}{\boldsymbol{\gamma}}^\text{T}{\boldsymbol{E}}_c^\text{T}{\boldsymbol{P}}{\boldsymbol{E}}_c{\boldsymbol{\gamma}}\;- \\ &\frac1\zeta\sum_{i=1}^3(h_i{\boldsymbol{N}}(\chi_i)-1)\dot{\chi}_i\leq\\ &-\alpha V-\frac1\zeta\sum_{i=1}^3(h_i{\boldsymbol{N}}(\chi_i)-1)\dot{\chi}_i+\rho\\[-1pt] \end{split} $$ (27) 式中, $ \alpha=-\lambda_{\max}({\boldsymbol{A}}_c+{\boldsymbol{P}}^{-1}{\boldsymbol{A}}_c^{\text{T}}{\boldsymbol{P}}+\varepsilon {\boldsymbol{I}}_6),\;\ \rho= \lambda_{\max}\times (({\boldsymbol{E}}_c^{\text{T}}{\boldsymbol{PE}}_c)\big/{2\varepsilon})\gamma_{\max}^2 $, $ \gamma_{\max} $是系统非线性估计误差$ \left\|\gamma\right\| $的上界; 因为LMI $ {\boldsymbol{PA}}_c+{\boldsymbol{A}}_c^\text{T}{\boldsymbol{P}}+\varepsilon {\boldsymbol{P}}<{\bf{0}},\;\ {\boldsymbol{P}}> {\bf{0}} $, 则$ {\boldsymbol{A}}_{c}+{\boldsymbol{P}}^{-1}{\boldsymbol{A}}_{c}^{\text{T}}{\boldsymbol{P}}+\varepsilon{\boldsymbol{I}}_{6} $的特征值为负实数, 故有$ \alpha>0 $.
接着, 基于定理2, 求解式(27)中的不等式, 其中, $ f=-1\big/{\zeta}\sum_{i=1}^3(h_i{\boldsymbol{N}}(\chi_i)-1)\dot{\chi}_i+\rho $.
$$ \begin{split}V(t)\leq\;&\text{e}^{-\alpha t}V(0)+\frac{\rho}{\alpha}(1-\text{e}^{-\alpha t})\;+\\&\text{e}^{-\alpha t}\int_0^t-\frac1\zeta\sum_{i=1}^3(h_i{\boldsymbol{N}}(\chi_i)-\text{l})\dot{\chi}_i\text{e}^{-\alpha\tau}\text{d}{{\tau}}\end{split} $$ (28) 显然, ${{\rm{e}}}^{-\alpha t}V(0)+\rho/{\alpha}(1-{{\rm{e}}}^{-\alpha t})$在区间$ [0,\;t_f] $上有界, 时变参数$ h_i,\;\ i=1,\;2,\;3 $取值范围在区间$ [h_{\min}, h_{\max}] $, $ \alpha>0 $, 根据定理1, 推知$\int_0^t-1/\zeta\sum_{i=1}^3 (h_i\times \boldsymbol{N}(\chi_i)-1)\dot{\chi}_i{\rm{e}}^{-\alpha\tau}\text{d}{\tau}\chi_i$和 $ V(t) $在区间 $ [0,\;t_f] $上均有界.
为确定系统收敛域, 设$ \mu_{\max} $, $ V_{\max} $分别为$\int_0^t- (1\big/{\zeta})\sum_{i=1}^3(h_i{\boldsymbol{N}}(\chi_i)-1)\dot{\chi}_i{{\rm{e}}}^{-\alpha\tau}\text{d}{{\tau}}$和 $ V(t) $的上确界, 则由式(28)可得:
$$ 0\leq V(t)\leq\kappa_1+\text{e}^{-\alpha t}\left(V(0)-\kappa_2\right)=V_{\max} $$ (29) 其中, $ \kappa_{1}=\kappa_{2}+\mu_{\max}\text{e}^{-\alpha t} $, $ \kappa_{2}=\rho/\alpha $. 这意味着姿态变量 ${\boldsymbol{X}}\;\;=\;\;[{\boldsymbol{\sigma}}^\text{T},\;\;{\boldsymbol{\dot{\sigma}}}^\text{T}]^\text{T}$ 的收敛域为 $ \left\|{\boldsymbol{X}}\right\|\;\leq \sqrt{2V_{\max}/\lambda_{\min}({\boldsymbol{P}})} $, 即当$ t\to\infty $系统状态最终收敛到一个平衡点附近的约束集$ {\boldsymbol{\Omega}}_{\boldsymbol{X}}\;=\;\{{\boldsymbol{X}}|\left\|{\boldsymbol{X}}\right\|\;\leq \sqrt{2\kappa_2/\lambda_{\min}\left({\boldsymbol{P}}\right)}\} $.
□ 4. 数学仿真及分析
本节将所设计的鲁棒饱和姿态控制方法应用于挠性航天器的姿态控制, 通过仿真验证该控制方法的有效性.
由图1所示航天器的材料特性与几何参数, 基于有限元方法进行动力学仿真, 前八阶模态振型如图3(a) ~ 3(h)所示. 其中, 一阶、二阶分别为$ xz $平面外绕$ y $轴的一阶横向弯曲(挥振)对称振型与$ yz $平面内绕$ x $轴的一阶横向弯曲(摆振)对称振型; 三阶、四阶分别为$ xz $平面外绕$ y $轴的一阶挥振反对称振型与$ yz$平面内绕$x$轴的一阶摆振反对称振型; 五阶、六阶分别为$ xz $平面外绕$ y $轴的二阶挥振的对称振型与反对称振型; 七阶、八阶分别为绕$ z$轴的一阶扭转反对称振型与对称振型.
考虑到低阶模态能量集度高, 基于模态截断法选取前四阶模态作为建立挠性航天器刚柔耦合动力学模型, 其合理性在后续仿真中将有所展示. 航天器所得动力学参数及外部环境干扰如表1所示.
表 1 仿真参数Table 1 Simulation parameters物理参数 值 转动惯量矩阵$ (\text{kg}\cdot\text{m}^{2}) $ ${\boldsymbol{J} }=\text{diag}\{40,\;150,\;160\}$ 耦合矩阵 $ {\boldsymbol{\delta}}= \left[\begin{array}{*{20}{r}} 1.352\ 3 & 1.278\ 4 & 2.155\ 3\\-1.151\ 9 & 1.017\ 6 & -1.272\ 4\\2.216\ 7 & 1.589\ 1 & -0.832\ 4\\1.236\ 4 & -1.653\ 7 & 1.225\ 1\end{array}\right] $ 挠性模态数 $ N=4 $ 固有频率(rad/s) $ {\boldsymbol{\Lambda}}=[1.20, \; 2.48 ,\; 3.37, \; 7.47] $ 阻尼比 $ \xi_1 = \xi_2 = \xi_3 = \xi_4 = 0.01 $ 外部干扰矩阵(N·m) $ d=10^{-4}\begin{bmatrix}3\cos(0.1t)+4\\1.5\sin(0.1t)+3\cos(0.1t)\\3\sin(0.1t)+1\end{bmatrix} $ 考虑执行机构的最大输出力矩为0.15 N·m, 欧拉角初值设置为[43.7, 38.5, 8.4] (°)、角速度初值为[0, 0, 0] (°/s), 模态位移速度初值[$ {\boldsymbol{\eta}} $, $ \dot{{\boldsymbol{\eta}}} $] = [$ {\bf{0}} $, $ {\bf{0}} $], 目标姿态、角速度、模态位移速度均为$ {\bf{0}} $. 本文提出的控制器对比PD + 干扰前馈补偿[2]、终端滑模控制(Terminal sliding mode control, TSMC)[5]、Zhao等[25]提出的综合全驱系统方法与扩展状态观测器(不对系统惯性耦合项进行观测)的姿态控制方法(Fully actuated system control based on extended state observer, FAESO)仿真结果如图4 ~ 图7所示, 各控制器参数设置如表2所示.
表 2 各控制器参数Table 2 Parameters of controllers控制方法 控制律 控制参数 PD + 前馈补偿 $ {\boldsymbol{u}} = {{\boldsymbol{B}}^{ - 1}}({{\boldsymbol{A}}_p}{\boldsymbol{\sigma}} {\rm{ + }}{{\boldsymbol{A}}_d}\dot{{\boldsymbol{\sigma}}} ) - {{\boldsymbol{B}}^{ - 1}}({{\boldsymbol{f}}_{\rm{2}}} + {{\boldsymbol{f}}_{\rm{3}}}) $ $\begin{array}{l} { {\boldsymbol{A} }_p} = {\rm{diag}}\{ - 0 .12 ,\; - 0.12 ,\; - 0.12 \}\\ { {\boldsymbol{A} }_d} = {\rm{diag\{ - 0.7 ,\; - 0.7,\; - 0.7} }\} \end{array}$ TSMC $ {\boldsymbol{u}} = - { {{\beta q} \over p}}{{\boldsymbol{B}}^{ - 1}}{\dot{{\boldsymbol{\sigma}}} ^{{\rm{(2}} - { {p \over q}}{\rm{)}}}} - {{\boldsymbol{B}}^{ - 1}}({{\boldsymbol{f}}_{\rm{2}}} + {{\boldsymbol{f}}_{\rm{3}}}) + \varepsilon \text{sgn} (s) $ $ \beta = 0.8,\;\varepsilon = 0.002,\; p = 5,\;q = 3 $ FAESO $ {\boldsymbol{u}} {\rm{ = }}{{\boldsymbol{B}}^{ - 1}}{\boldsymbol{Z}}{{\boldsymbol{F}}^2}{{\boldsymbol{V}}^{ - 1}}{[\sigma, {\rm{ }}\dot \sigma ]^{\rm{T}}} - {{\boldsymbol{B}}^{ - 1}}{\rm{(}}{\hat {{\boldsymbol{z}}}_3}{\rm{ + }}{{\boldsymbol{f}}_1}{\rm{)}} $ $\begin{array}{c} {\boldsymbol{F} } = \left[ \begin{array}{l} {\rm{diag} }\{ - 0.09,\; - 0.1,\; - 0.08\}\quad\ { {\boldsymbol{O} }_{3 \times 3} } \\ { {\boldsymbol{O} }_{3 \times 3{\rm{ } } } }\quad {\rm{diag} }\{ - 0.35,\; - 0.4,\; - 0.35\} \end{array} \right] \\ {\boldsymbol{Z} } = [{ {\boldsymbol{I} }_{3 \times 3} },\;\ { {\boldsymbol{I} }_{3 \times 3} }] \\ \left[\beta_1,\; \beta_2,\; \beta_3\right]=\left[50,\; 150,\; 500\right] \end{array}$ 本文方法 式(15)$ \sim $式(17) $\begin{array}{c} {\boldsymbol{F} } = \left[ \begin{array}{l} {\rm{diag} }\{ - 0.09,\; - 0.1,\; - 0.08\}\; { {\boldsymbol{O} }_{3 \times 3} }\\{ {\boldsymbol{O} }_{3 \times 3{\rm{ } } } }\; {\rm{diag} }\{ - 0.35,\; - 0.4,\; - 0.35\} \end{array} \right]\\{\boldsymbol{P} } = \left[ {\begin{array}{*{20}{c} } { {\rm{diag} }\{35,\;35,\;35\}\; {\rm{diag} }\{80,\;80,\;70\}}\\{ {\rm{diag} }\{80,\;80,\;80\}\;{\rm{diag} }\{150,\;150,\;150\}} \end{array} } \right]\\kv = 0.8,\;k = 0.6,\;\varsigma = 0.6,\; {\boldsymbol{Z} } = [{ {\boldsymbol{I} }_{3 \times 3} }{\rm{ } },\; { {\boldsymbol{I} }_{3 \times 3} }],\; {\boldsymbol{\chi} } (0){\rm{ = } }{\bf{0} }\\ \left[\beta_1,\; \beta_2,\; \beta_3\right]=\left[50,\; 150,\; 500\right],\;{ {\boldsymbol{z} }_{\rm{1} } } = {\bf{0} },\; { {\boldsymbol{z} }_{\rm{2} } } = {\bf{0} },\; { {\boldsymbol{z} }_{\rm{3} } } = {\bf{0} } \end{array}$ 图4 ~ 图5分别为欧拉角和姿态角速度的控制时间响应. 显然, 相较其他控制方法, 本文提出的控制器能够实现姿态的快速收敛, 且超调较小, 系统具有较好的动态特性. FAESO控制器次之. 在以上控制方法下的闭环姿态控制系统的三轴姿态精度均优于$ \left[0.015,\;1 \times {10^{-4}},\;7 \times {10^{-4}}\right]$ (°), 稳定度均优于$ \left[2 \times {10^{-3}},\;2 \times {10^{-5}},\;1.5 \times {10^{-4}}\right] $ (°/s).
图6为控制力矩$ {\boldsymbol{u}} $的时间响应. 由于惯性作用, 控制机构存在启动滞后性、姿态调整滞后性. 一方面, 该控制策略能有效减少由于执行器机构机械惯性引起的滞后性对系统性能的影响. 另一方面, 本文所提出的控制器考虑了饱和问题, 但相较FAESO算法, 控制输出集中于控制起始阶段, 可实现短时间内的控制效能最大化, 能够有效降低饱和约束所产生的影响, 明显缩短了饱和发生的时间. 另外, 在姿态机动过程中至机动末各控制方法下的控制力矩都存在抖振, 这很大程度上与存在的挠性振动与外部干扰有关, 但本方法产生的实际输出转矩相对平滑可行, 改善了系统的动态性能, 且能较快地完成控制任务.
图7为前四阶振动模态控制响应曲线. 显然, 四种控制器作用下挠性模态均能达到较小的稳态误差, 但以FAESO与本文方法稳态精度最高, 四阶模态精度均优于$ [6 \times {10^{-5}},\;2.5 \times {10^{-6}},\;2.0 \times {10^{-7}}, 1.5 \times 1{0^{-8}}] $, 在所提出控制器作用下的挠性模态具有较快时间响应, 振动衰减效果较为显著. 由此可见, 本文控制方法具有明显的挠性附件振动抑制效果. 另外, 相较其他控制器, 本控制器作用下前四阶模态响应在控制开始较大, 这与系统实际输出力矩有关, 控制开始执行机构满载运行, 与航天器本体姿态运动一起作用到挠性附件上, 引起较大的模态响应.
由图7可知, 一阶模态响应波动较大, 一阶模态对应平面$ xy $外绕$ y $轴的挥舞振动, 该方向与刚体偏航运动方向一致, 且能量集度较高的低阶模态振动最能表征物理坐标下的振动情况. 因此在刚柔耦合作用下, 挠性附件振动势必对航天器绕$ z $轴的姿态机动产生影响, 故一阶模态波动范围大的同时表现为在本文控制器作用下的系统绕$ z $轴的姿态相较其他方向存在较大波动, 如图4(c)所示. 此外, 随着模态阶数的增加, 模态响应呈跨量级式衰减, 说明取前四阶模态作为航天器动力学模型具有合理性.
挠性附件带来的挠性非线性对刚体姿态运动的扰动情况, 即刚柔耦合引起的非线性振动${\boldsymbol{f}}_2 $如图8所示, 起初挠性附件由于受航天器初始姿态作用对航天器的干扰力矩高达0.16 N·m, 在所提出控制器作用下, 挠性非线性逐渐衰减, 后续对航天器干扰较小.
努斯鲍姆函数的自适应参数$ {\boldsymbol{\chi}} $的时间响应如图9所示, 由式(17)可知, $ {\boldsymbol{\chi}} $与系统状态、控制律输出有关. 由图9所示, $ {\boldsymbol{\chi}} $的每个元素都是有界的, 并且单独收敛到某个具体常数, 即系统状态可控并趋于稳定状态, 控制律输出满足执行机构约束, 并趋于0.
5. 结束语
本文针对挠性航天器姿态机动与振动抑制问题, 提出了一种鲁棒饱和姿态控制方法. 该控制器的主要创新点为: 1)基于全驱系统理论采用直接参数法进行控制器设计, 将非线性姿态动力学转化为具有期望特征结构的线性定常系统, 有效缓解FAS理论的高维设计自由度带来的艰难选择. 2)利用努斯鲍姆增益技术处理输入饱和的影响, 使得系统实际控制力矩相对光滑可行. 3)引入非线性扩展观测器估计系统挠性振动、外部扰动等非线性项并进行有效补偿, 避免使用振动测量与控制所需应变片与作动器, 提高了系统可靠性并有效地节省星上空间资源. 考虑外界干扰、输入饱和、挠性振动及挠性模态不可测, 本文基于全驱系统理论, 重点研究存在外界干扰、输入饱和及挠性振动且挠性模态不可测等问题的带有大型挠性附件的航天器姿态镇定问题, 为后续研究关于挠性航天器的饱和姿态动态跟踪奠定理论基础, 为未来带有大型挠性附件航天器的姿态控制探索新方法.
-
表 1 仿真参数
Table 1 Simulation parameters
物理参数 值 转动惯量矩阵$ (\text{kg}\cdot\text{m}^{2}) $ ${\boldsymbol{J} }=\text{diag}\{40,\;150,\;160\}$ 耦合矩阵 $ {\boldsymbol{\delta}}= \left[\begin{array}{*{20}{r}} 1.352\ 3 & 1.278\ 4 & 2.155\ 3\\-1.151\ 9 & 1.017\ 6 & -1.272\ 4\\2.216\ 7 & 1.589\ 1 & -0.832\ 4\\1.236\ 4 & -1.653\ 7 & 1.225\ 1\end{array}\right] $ 挠性模态数 $ N=4 $ 固有频率(rad/s) $ {\boldsymbol{\Lambda}}=[1.20, \; 2.48 ,\; 3.37, \; 7.47] $ 阻尼比 $ \xi_1 = \xi_2 = \xi_3 = \xi_4 = 0.01 $ 外部干扰矩阵(N·m) $ d=10^{-4}\begin{bmatrix}3\cos(0.1t)+4\\1.5\sin(0.1t)+3\cos(0.1t)\\3\sin(0.1t)+1\end{bmatrix} $ 表 2 各控制器参数
Table 2 Parameters of controllers
控制方法 控制律 控制参数 PD + 前馈补偿 $ {\boldsymbol{u}} = {{\boldsymbol{B}}^{ - 1}}({{\boldsymbol{A}}_p}{\boldsymbol{\sigma}} {\rm{ + }}{{\boldsymbol{A}}_d}\dot{{\boldsymbol{\sigma}}} ) - {{\boldsymbol{B}}^{ - 1}}({{\boldsymbol{f}}_{\rm{2}}} + {{\boldsymbol{f}}_{\rm{3}}}) $ $\begin{array}{l} { {\boldsymbol{A} }_p} = {\rm{diag}}\{ - 0 .12 ,\; - 0.12 ,\; - 0.12 \}\\ { {\boldsymbol{A} }_d} = {\rm{diag\{ - 0.7 ,\; - 0.7,\; - 0.7} }\} \end{array}$ TSMC $ {\boldsymbol{u}} = - { {{\beta q} \over p}}{{\boldsymbol{B}}^{ - 1}}{\dot{{\boldsymbol{\sigma}}} ^{{\rm{(2}} - { {p \over q}}{\rm{)}}}} - {{\boldsymbol{B}}^{ - 1}}({{\boldsymbol{f}}_{\rm{2}}} + {{\boldsymbol{f}}_{\rm{3}}}) + \varepsilon \text{sgn} (s) $ $ \beta = 0.8,\;\varepsilon = 0.002,\; p = 5,\;q = 3 $ FAESO $ {\boldsymbol{u}} {\rm{ = }}{{\boldsymbol{B}}^{ - 1}}{\boldsymbol{Z}}{{\boldsymbol{F}}^2}{{\boldsymbol{V}}^{ - 1}}{[\sigma, {\rm{ }}\dot \sigma ]^{\rm{T}}} - {{\boldsymbol{B}}^{ - 1}}{\rm{(}}{\hat {{\boldsymbol{z}}}_3}{\rm{ + }}{{\boldsymbol{f}}_1}{\rm{)}} $ $\begin{array}{c} {\boldsymbol{F} } = \left[ \begin{array}{l} {\rm{diag} }\{ - 0.09,\; - 0.1,\; - 0.08\}\quad\ { {\boldsymbol{O} }_{3 \times 3} } \\ { {\boldsymbol{O} }_{3 \times 3{\rm{ } } } }\quad {\rm{diag} }\{ - 0.35,\; - 0.4,\; - 0.35\} \end{array} \right] \\ {\boldsymbol{Z} } = [{ {\boldsymbol{I} }_{3 \times 3} },\;\ { {\boldsymbol{I} }_{3 \times 3} }] \\ \left[\beta_1,\; \beta_2,\; \beta_3\right]=\left[50,\; 150,\; 500\right] \end{array}$ 本文方法 式(15)$ \sim $式(17) $\begin{array}{c} {\boldsymbol{F} } = \left[ \begin{array}{l} {\rm{diag} }\{ - 0.09,\; - 0.1,\; - 0.08\}\; { {\boldsymbol{O} }_{3 \times 3} }\\{ {\boldsymbol{O} }_{3 \times 3{\rm{ } } } }\; {\rm{diag} }\{ - 0.35,\; - 0.4,\; - 0.35\} \end{array} \right]\\{\boldsymbol{P} } = \left[ {\begin{array}{*{20}{c} } { {\rm{diag} }\{35,\;35,\;35\}\; {\rm{diag} }\{80,\;80,\;70\}}\\{ {\rm{diag} }\{80,\;80,\;80\}\;{\rm{diag} }\{150,\;150,\;150\}} \end{array} } \right]\\kv = 0.8,\;k = 0.6,\;\varsigma = 0.6,\; {\boldsymbol{Z} } = [{ {\boldsymbol{I} }_{3 \times 3} }{\rm{ } },\; { {\boldsymbol{I} }_{3 \times 3} }],\; {\boldsymbol{\chi} } (0){\rm{ = } }{\bf{0} }\\ \left[\beta_1,\; \beta_2,\; \beta_3\right]=\left[50,\; 150,\; 500\right],\;{ {\boldsymbol{z} }_{\rm{1} } } = {\bf{0} },\; { {\boldsymbol{z} }_{\rm{2} } } = {\bf{0} },\; { {\boldsymbol{z} }_{\rm{3} } } = {\bf{0} } \end{array}$ -
[1] 杨广华, 王强, 陈国玖, 秘倩. 美国“星链”低轨星座军事应用前景探析. 中国航天, 2022 (9): 60−63Yang Guang-Hua, Wang Qiang, Chen Guo-Jiu, Mi Qian. Future military application of the US starlink LEO constellation. Aerospace China, 2022 (9): 60−63 [2] Gu D K, Liu G P, Duan G R. Parametric control to a type of quasi-linear second-order systems via output feedback. International Journal of Control, 2019, 92(2): 291−302 doi: 10.1080/00207179.2017.1350885 [3] Yao Q J, Li Q, Huang J C, Jahanshahi H. PDE-based prescribed performance adaptive attitude and vibration control of flexible spacecraft. Aerospace Science and Technology, 2023, 141: Article No. 108504 doi: 10.1016/j.ast.2023.108504 [4] Li A, Liu M, Shi Y. Adaptive sliding mode attitude tracking control for flexible spacecraft systems based on the Takagi-Sugeno fuzzy modelling method. Acta Astronautica, 2020, 175: 570−581 doi: 10.1016/j.actaastro.2020.05.041 [5] Ghorbani H, Vatankhah R, Farid M. Adaptive nonsingular fast terminal sliding mode controller design for a smart flexible satellite in general planar motion. Aerospace Science and Technology, 2021, 119: Article No. 107100 doi: 10.1016/j.ast.2021.107100 [6] Hasan M N, Haris M, Qin S Y. Flexible spacecraft's active fault-tolerant and anti-unwinding attitude control with vibration suppression. Aerospace Science and Technology, 2022, 122: Article No. 107397 doi: 10.1016/j.ast.2022.107397 [7] Xiao Y, de Ruiter A, Ye D, Sun Z W. Adaptive fault-tolerant attitude tracking control for flexible spacecraft with guaranteed performance bounds. IEEE Transactions on Aerospace and Electronic Systems, 2022, 58(3): 1922−1940 doi: 10.1109/TAES.2021.3123295 [8] 张慧凤, 董乐伟, 魏新江. 一类带有输入饱和的随机系统抗干扰控制. 自动化学报, 2021, 47(6): 1453−1459Zhang Hui-Feng, Dong Le-Wei, Wei Xin-Jiang. Anti-disturbance control for a class of stochastic systems with input saturation. Acta Automatica Sinica, 2021, 47(6): 1453−1459 [9] 高阳, 吴文海, 王子健. 具有输入约束和输出噪声的不确定系统级联线性自抗扰控制. 自动化学报, 2022, 48(3): 843−852Gao Yang, Wu Wen-Hai, Wang Zi-Jian. Cascaded linear active disturbance rejection control for uncertain systems with input constraint and output noise. Acta Automatica Sinica, 2022, 48(3): 843−852 [10] 彭秀艳, 贾书丽, 张彪. 一类具有执行器饱和的非线性系统抗饱和方法研究. 自动化学报, 2016, 42(5): 798−804Peng Xiu-Yan, Jia Shu-Li, Zhang Biao. Anti-saturation method for a class of nonlinear systems with actuator saturation. Acta Automatica Sinica, 2016, 42(5): 798−804 [11] Piltan F, Mirzaei M, Shahriari F, Nazari I, Emamzadeh S. Design baseline computed torque controller. International Journal of Engineering, 2012, 6(3): 129−141 [12] Sancak K V, Bayraktaroglu Z Y. Nonlinear computed torque control of 6-dof parallel manipulators. International Journal of Control, Automation and Systems, 2022, 20(7): 2297−2311 doi: 10.1007/s12555-021-0198-6 [13] 段广仁. 高阶系统方法——II. 能控性与全驱性. 自动化学报, 2020, 46(8): 1571−1581Duan Guang-Ren. High-order system approaches: II. Controllability and full-actuation. Acta Automatica Sinica, 2020, 46(8): 1571−1581 [14] Duan G R. High-order fully actuated system approaches: Part III. Robust control and high-order backstepping. International Journal of Systems Science, 2021, 52(5): 952−971 doi: 10.1080/00207721.2020.1849863 [15] Duan G R. High-order fully actuated system approaches: Part IV. Adaptive control and high-order backstepping. International Journal of Systems Science, 2021, 52(5): 972−989 doi: 10.1080/00207721.2020.1849864 [16] Duan G R. High-order fully actuated system approaches: Part V. Robust adaptive control. International Journal of Systems Science, 2021, 52(10): 2129−2143 doi: 10.1080/00207721.2021.1879964 [17] Duan G R. High-order fully actuated system approaches: Part VIII. Optimal control with application in spacecraft attitude stabilisation. International Journal of Systems Science, 2022, 53(1): 54−73 doi: 10.1080/00207721.2021.1937750 [18] Xiao F Z, Chen L Q. Attitude control of spherical liquid-filled spacecraft based on high-order fully actuated system approaches. Journal of Systems Science and Complexity, 2022, 35(2): 471−480 doi: 10.1007/s11424-022-2055-y [19] Liu X Q, Chen M Y, Sheng L, Zhou D H. Adaptive fault-tolerant control for nonlinear high-order fully-actuated systems. Neurocomputing, 2022, 495: 75−85 doi: 10.1016/j.neucom.2022.04.129 [20] Hu Q L, Shao X D, Guo L. Adaptive fault-tolerant attitude tracking control of spacecraft with prescribed performance. IEEE/ASME Transactions on Mechatronics, 2018, 23(1): 331−341 doi: 10.1109/TMECH.2017.2775626 [21] Chen Z Y. Nussbaum functions in adaptive control with time-varying unknown control coefficients. Automatica, 2019, 102: 72−79 doi: 10.1016/j.automatica.2018.12.035 [22] Ioannou P A, Sun J. Robust Adaptive Control. Upper Saddle River: Prentice Hall, 1996. [23] Feng Z Y, Liu M, Cao X B. A fully-actuated system approach for spacecraft attitude control with input saturation. In: Proceedings of the 2nd Conference on Fully Actuated System Theory and Applications (CFASTA). Qingdao, China: IEEE, 2023. [24] 段广仁. 高阶系统方法——I. 全驱系统与参数化设计. 自动化学报, 2020, 46(7): 1333−1345Duan Guang-Ren. High-order system approaches: I. Fully-actuated systems and parametric designs. Acta Automatica Sinica, 2020, 46(7): 1333−1345 [25] Zhao Q, Duan G R. Fully actuated system approach for 6DOF spacecraft control based on extended state observer. Journal of Systems Science and Complexity, 2022, 35(2): 604−622 doi: 10.1007/s11424-022-1498-5 -