基于修正伪逆矩阵的冗余飞轮组合的角动量分配方法与流程

本发明涉及一种冗余飞轮组合的角动量分配方法，具体是指一种基于修正伪逆矩阵的冗余飞轮组合的角动量分配方法。

背景技术：

当前，基于角动量交换原理输出力矩的反作用飞轮是三轴姿态稳定卫星上使用最多的姿态控制执行机构。然而，现代航天任务对卫星的寿命提出了更高的要求，因此为了提高系统的可靠性，通常都会在卫星上安装多于3台的飞轮，构成冗余飞轮组合作为姿控系统的执行机构。由于飞轮数量大于输出力矩的维数，所以冗余飞轮组合在输出给定力矩的同时，还可以根据指令将飞轮转速偏置在最佳的工作转速附近。目前，冗余飞轮组合都是采用在伪逆分配算法中加入零运动的方式以改变飞轮的中心工作转速，例如国内专利《一种冗余飞轮角动量自主管理方法》就介绍了如何计算冗余飞轮组合的零运动以实现飞轮的转速偏置。

零运动的引入虽然可以改变飞轮组合的中心工作转速，但是由于每台飞轮角动量存在上限，所以这样做的结果势必会降低飞轮组合的角动量包络。以最简单的两台飞轮输出一维方向的力矩为例，假设单台飞轮输出力矩为0.2nm，角动量上限为25nms。当要求飞轮组合在正方向输出0.2nm的常值力矩时，如果初始状态下每台飞轮的角动量均为0nms，在采用伪逆加零运动的分配算法下，可以连续输出250s。但是，如果初始状态下，飞轮由于零运动的作用，已经分别偏置在了正负15nms的角动量下，则只能连续输出100s。上述两种情况下飞轮组合的初始合成角动量都是0nms，但在零运动偏置的情况下输出等值力矩的时间显著下降，表明零运动配合伪逆矩阵进行冗余飞轮角动量分配时，如果单台飞轮初始角动量不为零，会降低该飞轮角动量的利用率，导致飞轮组合在输出力矩方向的角动量仍有剩余的情况下，不能输出给定的力矩。

基于上述，亟需提出一种基于修正伪逆矩阵的冗余飞轮组合的角动量分配方法，可以根据冗余飞轮组合当前的状态，对伪逆矩阵进行实时的修正，以此抑制零运动对单台飞轮角动量利用率的影响，从而有效解决现有技术中存在的缺点和限制。

技术实现要素：

本发明的目的是提供一种基于修正伪逆矩阵的冗余飞轮组合的角动量分配方法，可以有效解决冗余飞轮组合在引入零运动后，单台飞轮角动量利用率显著下降的问题，可显著延长给定力矩的输出时间。

为实现上述目的，本发明提供一种基于修正伪逆矩阵的冗余飞轮组合的角动量分配方法，包含以下步骤：

s1、基于伪逆矩阵结合零运动的飞轮角动量分配算法，根据三轴指令力矩，计算冗余飞轮组合中的每台飞轮的参考驱动力矩；

s2、根据冗余飞轮组合中的每台飞轮的当前角动量，以及冗余飞轮组合中的每台飞轮的参考驱动力矩，对伪逆矩阵进行修正；包含以下步骤：

s21、根据每台飞轮的参考驱动力矩、每台飞轮的当前角动量以及每台飞轮的角动量上限，构建饱和变量，评价每台飞轮接近角动量饱和的程度；

s22、利用饱和变量计算性能指标修正系数，并利用饱和变量和性能指标修正系数得到修正后的伪逆矩阵；

s3、根据修正后的伪逆矩阵，计算冗余飞轮组合中的每台飞轮的驱动力矩。

所述的s1中，具体包含以下步骤：

s11、在一个控制周期ts内，对冗余飞轮组合中的各个飞轮的角动量增量δh进行分解，得到各个飞轮的输出力矩的角动量增量δht和各个飞轮的零运动的角动量增量δhz；

s12、通过伪逆矩阵ct(cct)-1，计算各个飞轮的输出力矩的角动量增量δht：

δht＝ct(cct)-1mts；

其中，c为飞轮安装矩阵；m为姿态控制器给出的当前控制周期需要施加的控制力矩；ts为姿态控制器的控制周期；

s13、计算各个飞轮的零运动的角动量增量δhz；

s14、计算各个飞轮的参考驱动力矩m*：

其中，分别为每台飞轮的参考驱动力矩。

所述的s21中，构建饱和变量αi的方法为：

其中，αi为第i台飞轮的饱和变量，代表第i台飞轮接近角动量饱和的程度；exp(…)为指数函数；min(…)为求最小值函数；sign(…)为符号函数；himax表示第i台飞轮的最大角动量；表示第i台飞轮的参考驱动力矩；hmi表示第i台飞轮的当前角动量；n为冗余飞轮组合中的飞轮数量。

所述的s22中，具体包含以下步骤：

s221、根据饱和变量αi计算性能指标修正系数η：

s222、根据饱和变量αi和性能指标修正系数η，构建性能指标加权系数矩阵r：

s223、根据性能指标加权系数矩阵r，对伪逆矩阵ct(cct)-1进行修正，得到修正后的伪逆矩阵rct(crct)-1。

所述的s3中，具体包含以下步骤：

s31、根据修正后的伪逆矩阵rct(crct)-1，计算修正后的各个飞轮的输出力矩的角动量增量

s32、根据修正后的各个飞轮的输出力矩的角动量增量基于伪逆矩阵结合零运动的飞轮角动量分配算法，计算各个飞轮的驱动力矩m：

其中，m＝[m1m2...mn]t，且m1，m2，...mn分别为每台飞轮的驱动力矩。

综上所述，本发明提供的基于修正伪逆矩阵的冗余飞轮组合的角动量分配方法，与现有技术相比，具有以下优点和有益效果：

1、本发明能够有效解决冗余飞轮组合在引入零运动后，单台飞轮角动量利用率显著下降的问题，可显著延长给定力矩的输出时间；

2、本发明无需反复迭代寻优，计算量小，占用星上资源少，易于实际应用；

3、本发明以解析的形式给出，可以实时进行计算，不存在数值不稳定的问题。

附图说明

图1为本发明中未使用修正伪逆矩阵的冗余飞轮组合的输出力矩的示意图；

图2为本发明中使用修正伪逆矩阵后的冗余飞轮组合的输出力矩的示意图；

图3为本发明中未使用修正伪逆矩阵的冗余飞轮组合的每台飞轮角动量的示意图；

图4为本发明中使用修正伪逆矩阵后的冗余飞轮组合的每台飞轮角动量的示意图；

图5为本发明中基于修正伪逆矩阵的冗余飞轮组合的角动量分配方法的流程图。

具体实施方式

以下结合图1～图5，通过优选实施例对本发明的技术内容、构造特征、所达成目的及功效予以详细说明。

首先说明的是，本发明中统一使用大写加粗向量作为冗余飞轮组合的合成角动量或输出/驱动力矩，统一使用小写加粗向量作为冗余飞轮组合内的每台单台飞轮的角动量或输出/驱动力矩。

如图5所示，为本发明提供的基于修正伪逆矩阵的冗余飞轮组合的角动量分配方法，包含以下步骤：

s1、按照现有技术中的基于伪逆矩阵结合零运动的飞轮角动量分配算法，根据三轴指令力矩，计算冗余飞轮组合中的每台飞轮的参考驱动力矩；

s2、根据冗余飞轮组合中的每台飞轮的当前角动量，以及s1中计算得到的冗余飞轮组合中的每台飞轮的参考驱动力矩，对伪逆矩阵进行修正，从而提高每台飞轮的角动量利用率；

s3、根据修正后的伪逆矩阵，计算冗余飞轮组合中的每台飞轮的驱动力矩。

所述的s1中，具体包含以下步骤：

s11、对冗余飞轮组合的合成角动量增量进行分解；

在一个控制周期ts内，冗余飞轮组合的合成角动量增量是由两部分组成，即：

δh＝δht+δhz

其中，表示冗余飞轮组合的输出力矩的角动量增量，表示冗余飞轮组合的零运动的角动量增量；表示多维实数矩阵空间，其中右上角的数字3×1则表示该矩阵的行数和列数(下文中的都是表示多维实数矩阵空间)；

在一个控制周期ts内，冗余飞轮组合中的各个飞轮的角动量增量δh是由两部分组成，即：

cδh＝cδht+cδhz

其中，c表示飞轮安装矩阵；δht表示各个飞轮的输出力矩的角动量增量；δhz表示各个飞轮的零运动的角动量增量；

需要说明的是，δh是以每台飞轮的角动量增量为元素组成的n维向量；δht是以每台飞轮的输出力矩的角动量增量为元素组成的n维向量；δhz是以每台飞轮的零运动的角动量增量为元素组成的n维向量；n表示冗余飞轮组合中的飞轮数量；

s12、计算各个飞轮的输出力矩的角动量增量δht；

根据冗余飞轮组合的输出力矩的角动量增量δht，通过伪逆矩阵ct(cct)-1，可计算得到：

δht＝ct(cct)-1δht＝ct(cct)-1mts

其中，m为姿态控制器给出的当前控制周期需要施加的控制力矩；ts为姿态控制器的控制周期；

s13、计算各个飞轮的零运动的角动量增量δhz；

设各个飞轮的中心角动量为hz，同样的，hz是以每台飞轮的中心角动量为元素组成的n维向量；令：

hz＝be

ce＝0

其中，b为人为设定的标量，选择依赖工程经验，设计原则为根据卫星所受干扰力矩在一个轨道周期内的变化和积累情况，尽量保证每台飞轮都工作在理想转速附近；e为飞轮安装矩阵c的零空间基解，设计的原则为e中各元素的大小尽量接近；

各个飞轮的零运动的角动量增量δhz可通过采用基于qr分解的方法求解如下一个有等式约束的二次最优问题得到：

s.t.cδhz＝0

其中，hm＝[hm1hm2...hmn]t为各个飞轮的当前角动量，hm1，hm2，...hmn分别为每台飞轮的当前角动量；

求解得到：

其中，k是控制增益矩阵；q2可通过以下方法求得；

对飞轮安装矩阵c进行qr分解：

其中，q为分解得到的酉矩阵，r为上三角矩阵；而酉矩阵q可以写成如下的分块矩阵形式：

其中，q1和q2均是列数为n/2的矩阵；

假设令控制增益矩阵k采用等比例限幅的形式表示，即

其中，δhmax是单台飞轮在一个控制周期内的角动量零运动的上限；k1,...,kn是零运动控制增益；diag(…)表示以括号内元素构成的对角线矩阵；max(…)是求最大值函数；

此处进行限幅处理是考虑到单台飞轮输出力矩存在上限，防止单台飞轮因零运动太大而导致该飞轮无法正常输出力矩。

s14、根据各个飞轮的输出力矩的角动量增量δht以及各个飞轮的零运动的角动量增量δhz，计算各个飞轮的参考驱动力矩；具体为：

m*＝δht/ts+δhz/ts

其中，m*为各个飞轮的参考驱动力矩，且m*是由每台飞轮的参考驱动力矩为元素组成的向量；分别为每台飞轮的参考驱动力矩。

所述的s2为本发明的核心步骤，该s2中具体包含以下步骤：

s21、以每台飞轮的参考驱动力矩为基础，结合每台飞轮的当前角动量以及每台飞轮的角动量上限，构建饱和变量α以评价每台飞轮接近角动量饱和的程度；

设各个飞轮的最大角动量为hmax，且hmax是由每台飞轮的最大角动量为元素组成的向量，则有：

hmax＝[h1maxh2max...hnmax]t

其中，h1max，h2max，...hnmax分别为每台飞轮的最大角动量；

根据第i台飞轮的参考驱动力矩和当前角动量hmi，取饱和变量αi代表第i台飞轮距离角动量饱和的接近程度，αi越大，代表第i台飞轮越接近角动量饱和；其中，i＝1,...,n，且n为冗余飞轮组合中的飞轮数量；

采用下式计算饱和变量αi：

其中，exp(…)为指数函数；min(…)为求最小值函数；sign(…)为符号函数；himax表示第i台飞轮的最大角动量；

s22、利用饱和变量α计算性能指标修正系数η，并利用饱和变量α和性能指标修正系数η得到修正后的伪逆矩阵；

求取各个飞轮的输出力矩的角动量增量δht的伪逆矩阵是通过求取如下最优问题得到的：

s.t.cδht＝mts

根据饱和变量α和性能指标修正系数η，将上述最优问题修正为：

s.t.cδht＝mts

其中，r为性能指标加权系数矩阵；且有：

根据性能指标加权系数矩阵r，对伪逆矩阵ct(cct)-1进行修正，得到修正后的伪逆矩阵rct(crct)-1。

所述的s3中，具体包含以下步骤：

s31、根据修正后的伪逆矩阵，计算修正后的各个飞轮的输出力矩的角动量增量

其中，m为姿态控制器给出的当前控制周期需要施加的控制力矩；ts为姿态控制器的控制周期；

s32、根据修正后的各个飞轮的输出力矩的角动量增量基于伪逆矩阵结合零运动的飞轮角动量分配算法，计算各个飞轮的驱动力矩m：

其中，m＝[m1m2...mn]t，且m1，m2，...mn分别为每台飞轮的驱动力矩。

本发明提供的基于修正伪逆矩阵的冗余飞轮组合的角动量分配方法，其有效性可通过仿真进行验证。如图1～图4所示，选取4台飞轮构成冗余飞轮组合，单台飞轮角动量上限为25nms。首先通过零运动将每台飞轮的角动量偏置在7.5nms附近，在500s时开始持续在x方向输出0.1nm的力矩。可以看到对于未使用修正伪逆矩阵的角动量分配方法，在690s时就有一台飞轮出现了角动量饱和的情况。而对于使用了修正伪逆矩阵的角动量分配方法，一直要到760s时才会出现单台飞轮的角动量饱和情况。因此，显而易见，本发明方法可以显著延长给定力矩的输出时间。

综上所述，本发明提供的基于修正伪逆矩阵的冗余飞轮组合的角动量分配方法，能够有效解决冗余飞轮组合在引入零运动后，单台飞轮角动量利用率显著下降的问题，同时可增加使用零运动后冗余飞轮组合角动量包络的范围。

尽管本发明的内容已经通过上述优选实施例作了详细介绍，但应当认识到上述的描述不应被认为是对本发明的限制。在本领域技术人员阅读了上述内容后，对于本发明的多种修改和替代都将是显而易见的。因此，本发明的保护范围应由所附的权利要求来限定。