CN103543750A - 基于二次型规划的在轨服务航天器推力分配优化方法 - Google Patents

Abstract

本发明提供了一种基于二次型规划的在轨服务航天器推力分配优化方法，采用基于二次型规划的二次分步优化方法，首先应用伪逆法寻找满足轨道与姿态控制要求的初始解，再使用二次型规划可以将求得的初始解修正到推力器能够提供的推力范围内。这种方法避免了单独采用二次型规划方法寻求初始解的盲目性和不精确性，使得期望控制量能够在执行机构间进行更好的分配，进而提高运算效率。

Description

基于二次型规划的在轨服务航天器推力分配优化方法

技术领域

本发明涉及一种航天器推力分配方法。 

背景技术

在现有的在轨服务任务中，航天器的推力分配，大都是根据推力器布局，预先制定分配列表。该方法的主要缺点是需要预先制定推力器分配列表，包括推力器故障时的分配列表，它需要占用大量的星上存储空间，且采用这种分配方式也无法实时地调整分配策略应对不可遇见的推力器故障。因此，有必要基于现有的优化理论，研究新型的推力分配优化算法，实现在空间操作任务中，对装配有多推力器系统的航天器进行有效，实时的推力分配，同时减小在控制过程中的燃料消耗。 

在空间操作任务中，执行在轨服务的主动航天器需同时控制主动航天器和目标航天器之间的相对位置和相对姿态。在这种情况下，控制执行机构由十几台到数十台推力器组成，之所以需要如此多的推力器，是因为在轨服务需要进行三轴相对位置控制和三轴相对姿态控制，且对控制精度的要求极高。通过公用一套推进系统实现轨道和姿态的一体化控制，则既能提高航天器的执行精度和对推力器故障的容错能力，又能节省部分硬件、节约燃料。 

为了确保航天器的在轨可靠运行，现有设计理念通常采用冗余配置系统，这使得由控制算法给出的期望控制量到推力器控制指令的分配方案并不唯一。控制分配方法是由控制律给出的期望控制量出发，在各类约束条件和最优目标下，将期望控制量在冗余配置的执行机构间进行分配，使执行机构实际控制输出尽可能与期望控制量相吻合的一种控制设计技术。控制分配算法的实质是一个约束条件下期望控制量到各个执行机构控制指令的非线性映射过程。基于优化的控制分配算法将控制分配问题转化为包括代价函数、等式约束和不等式约束的约束优化数学模型，然后利用各种优化算法对其进行求解。 

总的来说，当前还未能提出一种实时有效的航天器推力分配优化算法。 

发明内容

为了克服现有技术的不足，本发明提供一种基于二次型规划的在轨服务航天器推力分配优化方法，以实现在轨服务航天器在空间操作任务中的自主推力分配。 

本发明解决其技术问题所采用的技术方案包括以下步骤： 

步骤一、定义航天器本体坐标系Ox_by_bz_b的三个坐标轴与航天器三个惯性主轴重合，航天器上安装的推力器个数为n个，在航天器体坐标系Ox_by_bz_b下，各个推力器相对航天器质心的位置矢量矩阵为[d_1b,d_2b,…,d_nb]，其中d_ib=[x_i  y_i  z_i]^T表示d_ib在体坐标系Ox_by_bz_b三个轴上的分量大小；定义e_b=[e_x  e_y  e_z]^T为体坐标系的三个基矢，推力器所产生的单位推力矩阵为[e_1b,e_2b,…,e_nb]，其中e_ib=[cosα_icosβ_icosα_icosβ_isinα_i]^T为第i个推力器产生的单位推力在体坐标系三个轴上的分量；第i个推力器产生的推力大小为F_i，i＝1,2,…,n，且

为第i个推力器在一个执行周期内全部开机时所能产生的最大推力，则第i个推力器对航天器质心产生的作用力为： 

$U_i=e_b^T{U}_{\mathrm{ib}}=e_b^T{F}_i{e}_{\mathrm{ib}}---\left(1\right)$

故 

U_ib＝F_ie_ib    (2) 

其产生绕质心的作用力矩为： 

$T_i=e_b^T{T}_{\mathrm{ib}}=e_b^T{d}_{\mathrm{ib}}\×U_{\mathrm{ib}}---\left(3\right)$

故 

T_i＝(d_ib×e_ib)F_i                                   (4) 

设所有推力器产生的推力组成列阵F＝[F₁,F₂,…,F_n]^T，则其在航天器质心处合成的作用力矩可以表示为： 

$T=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{T}_i=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}(d_{\mathrm{ib}}\×e_{\mathrm{ib}})F_i=\mathrm{AF}---\left(5\right)$

所产生的作用力可以表示为： 

$U=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{U}_i=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{F}_i{e}_{\mathrm{ib}}=\mathrm{BF}---\left(6\right)$

式中，A为所有推力器的单位推力矢量对航天器的力矩矩阵： 

A＝[d_1b×e_1b,d_2b×e_2b,…,d_nb×e_nb]                    (7)B为所有推力器的单位推力矢量对航天器的力矩阵： 

B＝[e_1b,e_2b,…,e_nb]                                  (8) 

令C＝[T_c，U_c]^T _,D＝[A，B]^T，其中T_c为期望控制力矩，U_c为期望控制力，则推力分配的数学描述为： 

C=DF    (9) 

步骤二、考虑泛函： 

$J\left(F\right)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{F_i}^2=F^{T}F---\left(10\right)$

可知，推力F须满足的约束条件为： 

C=DF                 (11)其中，C＝[T_c;U_c],D＝[A;B]，T_c为期望控制力矩，U_c为期望控制力； 

得到推力F的表达式为： 

F＝D^T(DD^T)^-1C＝D⁺C         (12) 

其中，D⁺＝D^T(DD^T)^-1即为D的伪逆； 

对推力F进行修正，令 

F＝D⁺C+w          (13) 

其中w为修正变量，满足齐次线性方程： 

Dw＝0         (14) 

式(14)的解表示为： 

w＝k₁ξ₁+k₂ξ₂+…+k_n-6ξ_n-6＝ξk        (15) 

其中，ξ＝[ξ₁,ξ₂,…,ξ_n-6]为方程组的一个基础解系；k＝[k₁,k₂,…,k_n-6]^T，k_i为任意实数1≤i≤n-6,i∈N； 

考虑泛函： 

$\min J\left(F\right)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{F_i}^2=F^{T}F$

$={\left(\right[D^{+}C]+w)}^T\left(\right[D^{+}C]+w)---\left(16\right)$

$=\left\{{\left[D^{+}C\right]}^T\right[D^{+}C\left]\right\}+\{{\left[D^{+}C\right]}^{T}w+w^T[D^{+}C]+w^{T}w\}$

上述泛函等价于： 

$\min J\left(w\right)=\{{\left[D^{+}C\right]}^{T}w+w^T[D^{+}C]+w^{T}w\}=\frac{1}{2}{k}^T\mathrm{Hk}+f^{T}k---\left(17\right)$

其中f＝2(ξ^T[D⁺C])，H＝2ξ^Tξ； 再考虑F的有界性，亦即单个推力器推力范围限制：

令  $F^{\mathrm{\μ}}={[F_1^{\mathrm{\μ}},F_2^{\mathrm{\μ}},...,F_n^{\mathrm{\μ}}]}^T,$ 考虑式(13)则得： 

$\begin{array}{c}w=\mathrm{\ξk}\≤F^{\mathrm{\μ}}-D^{+}C\\ -w=-\mathrm{\ξk}\≤D^{+}C\end{array}---\left(18\right)$

令G＝[ξ;-ξ]^T,S＝[F^μ-D⁺C;D⁺C]^T，则约束条件(18)表示为：Gk≤S                                     (19) 

式(17)与式(19)构成二次型规划问题的一般形式，对该轨迹规划问题进行求解，从而得到最优解k^*，亦即得到满足轨道与姿态期望控制要求，且满足推力器约束的推力F^*，最终完成推力分配过程。 

本发明的有益效果是：由于采用基于二次型规划的二次分步优化方法。首先应用伪逆法寻找满足轨道与姿态控制要求的初始解，再使用二次型规划可以将求得的初始解修正到推力器能够提供的推力范围内。这种方法避免了单独采用二次型规划方法寻求初始解的盲目性和不精确性，使得期望控制量能够在执行机构间进行更好的分配，进而提高运算效率。 

附图说明

图1是第i个单位推力矢量在体系下的分量对位置变化示意图。 

具体实施方式

下面结合附图和实施例对本发明进一步说明，本发明包括但不仅限于下述实施例。 

本发明针对航天器上装有一套多台推力器组成的推进系统以期实现轨道和姿态一体化控制的情况。 

本发明由以下步骤实现： 

步骤一、建立控制分配问题的数学模型 

定义航天器本体坐标系Ox_by_bz_b的三个坐标轴与航天器三个惯性主轴重合，航天器上安装的推力器个数为n个。在航天器体坐标系Ox_by_bz_b下，它们相对航天器质心的位置矢量矩阵为[d_1b,d_2b,…,d_nb]，其中d_ib=[x_i  y_i  z_i]^T表示d_ib在体坐标系Ox_by_bz_b三个轴上的分量大小。定义e_b=[e_x  e_y  e_z]^T为体坐标系的三个基矢，推力器所产生的单位推力矩阵为[e_1b,e_2b,…,e_nb]，其中e_ib=[cosα_icosβ_icosα_icosβ_isinα_i]^T为第i个推力 器产生的单位推力在体坐标系三个轴上的分量。如图1所示。 

若第i个推力器才所产生的推力大小为F_i(i＝1,2,…,n)，且

为第i个推力器在一个执行周期内全部开机时所能产生的最大推力，则它对航天器质心产生的作用力为： 

$U_i=e_b^T{U}_{\mathrm{ib}}={e_b^{T}F}_i{e}_{\mathrm{ib}}---\left(1\right)$

故 

U_ib＝F_ie_ib    (2) 

其产生绕质心的作用力矩为： 

$T_i=e_b^T{T}_{\mathrm{ib}}=e_b^T{d}_{\mathrm{ib}}\×U_{\mathrm{ib}}---\left(3\right)$

故 

T_i＝(d_ib×e_ib)F_i           (4) 

设所有推力器产生的推力组成列阵F＝[F₁,F₂,…,F_n]^T，则其在航天器质心处合成的作用力矩可以表示为： 

$T=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{T}_i=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}(d_{\mathrm{ib}}\×e_{\mathrm{ib}})F_i=\mathrm{AF}---\left(5\right)$

所产生的作用力可以表示为： 

$U=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{U}_i=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{F}_i{e}_{\mathrm{ib}}=\mathrm{BF}---\left(6\right)$

式中，A为所有推力器的单位推力矢量对航天器的力矩矩阵： 

A＝[d_1b×e_1b,d_2b×e_2b,…,d_nb×e_nb]       (7)B为所有推力器的单位推力矢量对航天器的力矩阵： 

B＝[e_1b,e_2b,…,e_nb]       (8) 

令C＝[T_c，U_c]^T,D＝[A，B]^T，其中T_c为期望控制力矩，U_c为期望控制力。则推力分配的数学描述为： 

C=DF    (9) 

步骤二、基于二次型规划的推力分配优化算法 

考虑泛函： 

$J\left(F\right)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{F_i}^2=F^{T}F---\left(10\right)$

上式表示的泛函指标满足其所有元素的平方和最小。由式(9)可知，推力F须满足的约束条件为： 

C=DF        (11)其中，C＝[T_c;U_c],D＝[A;B]。T_c为期望控制力矩，U_c为期望控制力。 

根据矩阵理论，通过伪逆法得到推力F的表达式为： 

F＝D^T(DD^T)^-1C＝D⁺C       (12)其中，D⁺＝D^T(DD^T)^-1即为D的伪逆。 

利用伪逆法求取的推力能够满足航天器轨道与姿态控制要求，但它存在一定的缺陷，亦即不能保证求得的推力在推力器所能提供的推力大小范围内，这就需要对其进行修正。由于采用伪逆方法求得的初始解不一定满足推力等约束条件，故采用二次型规划方法将首次优化结果修正到推力器所能提供的推力范围内。 

令 

F＝D⁺C+w        (13)其中w为修正变量，满足齐次线性方程： 

Dw＝0         (14) 

将式(13)两端同时左乘D，在式(11)表示的约束条件下，可以得到式(14)。可见在不影响航天器轨道与姿态控制要求的前提下，适当选取w，可以将首次优化的F调节到推力器所能提供的推力范围内。式(14)的解可以表示为： 

w＝k₁ξ₁+k₂ξ₂+…+k_n-6ξ_n-6＝ξk        (15)其中，ξ＝[ξ₁,ξ₂,…,ξ_n-6]为方程组的一个基础解系；k＝[k₁,k₂,…,k_n-6]^T，k_i为任意实数(1≤i≤n-6,i∈N)。以下的任务是选取适当的系数k。 

考虑泛函： 

$\min J\left(F\right)=\underset{i=1}{\overset{n}{\mathrm{\Σ}}}{F_i}^2=F^{T}F$

$={\left(\right[D^{+}C]+w)}^T\left(\right[D^{+}C]+w)---\left(16\right)$

$=\left\{{\left[D^{+}C\right]}^T\right[D^{+}C\left]\right\}+\{{\left[D^{+}C\right]}^{T}w+w^T[D^{+}C]+w^{T}w\}$

由于上式第一项为定值，故上述泛函等价于： 

$\min J\left(w\right)=\{{\left[D^{+}C\right]}^{T}w+w^T[D^{+}C]+w^{T}w\}=\frac{1}{2}{k}^T\mathrm{HK}+f^{T}k---\left(17\right)$

其中f＝2(ξ^T[D⁺C])，H＝2ξ^Tξ。 

再考虑F的有界性，亦即单个推力器推力范围限制：令  $F^{\mathrm{\μ}}={[F_1^{\mathrm{\μ}},F_2^{\mathrm{\μ}},...,F_n^{\mathrm{\μ}}]}^T,$ 考虑式(13)则得： 

$\begin{array}{c}w=\mathrm{\ξk}\≤F^{\mathrm{\μ}}-D^{+}C\\ -w=-\mathrm{\ξk}\≤D^{+}C\end{array}---\left(18\right)$

令G＝[ξ;-ξ]^T,S＝[F^μ-D⁺C;D⁺C]^T，则约束条件(18)可以表示为： 

Gk≤S                                     (19) 

式(17)与式(19)构成二次型规划问题的一般形式，至此已将自主交会的轨迹规划问题转换为一个标准的二次型规划问题，可使用Matlab软件当中的quadprog函数对该轨迹规划问题进行求解，从而得到最优解k^*，亦即得到满足轨道与姿态期望控制要求，且满足推力器约束的推力F^*，最终完成推力分配过程。 

本发明方法的实例验证： 

1）单个推力器最大推力

推力器距质心的力臂为L＝0.4m； 

2）所有推力器的单位推力所产生的力矩矩阵为： 

$A=\left[\begin{array}{cccccccccccccccc}0& 0& 0& 0& 0& 0& 0& 0& -0.4& 0.4& 0& 0& 0.4& -0.4& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 0& 0& -0.4& 0.4& 0& 0& 0.4& -0.4\\ 0.4& -0.4& 0.4& -0.4& -0.4& 0.4& -0.4& 0.4& 0& 0& 0& 0& 0& 0& 0& 0\end{array}\right]$

3）而所有推力器的单位推力产生的力矩阵为： 

$B=\left[\begin{array}{cccccccccccccccc}1& 1& 0& 0& -1& -1& 0& 0& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 1& 1& 0& 0& -1& -1& 0& 0& 0& 0& 0& 0& 0& 0\\ 0& 0& 0& 0& 0& 0& 0& 0& 1& 1& 1& 1& -1& -1& -1& -1\end{array}\right];$

4）期望控制量列表为； 

表1期望控制量 

序号	期望控制量
		1	[0 0 0 0.01 0.01 0.02]T
2	[0.003 0.004 0.005 0 0 0]T
		3	[0.003 0.004 0.005 0.01 0.02 0.03]T
4	[0.004 0.001 0.003 0.01 0.02 0.01]T
		5	[0.004 0.001 0.012 0.08 0.02 0.01]T

 在上述给定条件下，采用伪逆法和二次型规划法结合的动态分配方法，对期望控制量进行分配，并与传统固定分配方案模式相比较，结果如下表所示： 

表2推力分配结果与燃料消耗对比 

从上表可以看出： 

1)本发明提出的分配算法能够在推力器配置的负载范围内有效完成航天器一体化控制任务； 

2)对于只有期望控制力或期望控制力矩情况1和2，算法给出的控制量与传统控制模式控制量的燃料消耗相同，而在其它情况下，本发明所设计的分配算法都能显著降低燃料消耗。 

Claims (1)

1.一种基于二次型规划的在轨服务航天器推力分配优化方法，其特征在于包括下述步骤： 

步骤一、定义航天器本体坐标系Ox_by_bz_b的三个坐标轴与航天器三个惯性主轴重合，航天器上安装的推力器个数为n个，在航天器体坐标系Ox_by_bz_b下，各个推力器相对航天器质心的位置矢量矩阵为[d_1b,d_2b,…,d_nb]，其中d_ib=[x_i y_i z_i]^T表示d_ib在体坐标系Ox_by_bz_b三个轴上的分量大小；定义e_b=[e_x e_y e_z]^T为体坐标系的三个基矢，推力器所产生的单位推力矩阵为[e_1b,e_2b,…，e_nb]，其中e_ib=[cosα_icosβ_i cosα_icosβ_i sinα_i]^T为第i个推力器产生的单位推力在体坐标系三个轴上的分量；第i个推力器产生的推力大小为F_i，i＝1,2,…,n，且

为第i个推力器在一个执行周期内全部开机时所能产生的最大推力，则第i个推力器对航天器质心产生的作用力为： 

故 

U_ib＝F_ie_ib    (2) 

其产生绕质心的作用力矩为： 

故 

T_i＝(d_ib×e_ib)F_i    (4) 

设所有推力器产生的推力组成列阵F＝[F₁,F₂,…,F_n]^T，则其在航天器质心处合成的作用力矩可以表示为： 

所产生的作用力可以表示为： 

式中，A为所有推力器的单位推力矢量对航天器的力矩矩阵： 

A＝[d_1b×e_1b,d_2b×e_2b,…，d_nb×e_nb]    (7) 

B为所有推力器的单位推力矢量对航天器的力矩阵： 

B＝[e_1b,e_2b,…,e_nb]    (8) 

令C＝[T_c，U_c]^T,D＝[A，B]^T，其中T_c为期望控制力矩，U_c为期望控制力，则推力分配的数学描述为： 

C=DF    (9) 

步骤二、考虑泛函： 

可知，推力F须满足的约束条件为： 

C=DF    (11) 

其中，C＝[T_c;U_c],D＝[A;B]，T_c为期望控制力矩，U_c为期望控制力； 

得到推力F的表达式为： 

F＝D^T(DD^T)^-1C＝D⁺C    (12) 

其中，D⁺＝D^T(DD^T)^-1即为D的伪逆； 

对推力F进行修正，令 

F＝D⁺C+w    (13) 

其中w为修正变量，满足齐次线性方程： 

Dw＝0    (14) 

式(14)的解表示为： 

w＝k₁ξ₁+k₂ξ₂+…+k_n-6ξ_n-6＝ξk    (15) 

其中，ξ＝[ξ₁,ξ₂,…,ξ_n-6]为方程组的一个基础解系；k＝[k₁,k₂,…,k_n-6]^T，k_i为任意实数1≤i≤n-6,i∈N； 

考虑泛函： 

上述泛函等价于： 

其中f＝2(ξ^T[D⁺C])，H＝2ξ^Tξ； 

再考虑F的有界性，亦即单个推力器推力范围限制：

令 

考虑式(13)则得： 

令G＝[ξ;-ξ]^T,S＝[F^μ-D⁺C;D⁺C]^T，则约束条件(18)表示为： 

Gk≤S    (19) 

式(17)与式(19)构成二次型规划问题的一般形式，对该轨迹规划问题进行求解，从而得到最优解k^*，亦即得到满足轨道与姿态期望控制要求，且满足推力器约束的推力F^*，最终完成推力分配过程。 

Images