CN102354218B - 一种航天器相对运动的采样控制方法 - Google Patents

Abstract

一种航天器相对运动的采样控制方法，它涉及一种航天器的采样控制方法。本发明为解决采用现有的航天器相对运动的采样控制方法忽略了数字控制器的处理周期和偏差，影响航天器轨道的精确性和安全性的问题。步骤A、建立航天器相对运动动力学模型；步骤B、对两个航天器相对状态进行采样；步骤C、利用步骤B中所述的扇形区域的上下边界线构造M和N矩阵；步骤D、求得相应的状态反馈控制律；步骤E、引入两个正定对称矩阵P和Q并定义如下李亚普诺夫泛函；步骤F、求得交会过程完成并且推力满足公式(3)上界约束条件；步骤G、利用MATLAB软件中线性矩阵不等式(LMI)工具箱求可行解。本发明的采样控制方法用于设计航天器控制器。

Description

一种航天器相对运动的采样控制方法

技术领域

本发明涉及一种航天器的采样控制方法。

背景技术

连续推力控制是一种重要的航天器轨道机动控制形式，在航天器自主交会、编队飞行、空间站停靠等多种航天器相对运动任务中获得广泛应用。

目前很多利用连续推力形式的轨道机动控制方法完全基于连续系统模型和连续控制器形式。但随着计算机技术的飞速发展，实际工程中采用的控制器多为数字信号形式的计算机系统。在这种系统中，控制过程需要利用采样器对航天器相对运动状态进行固定时间间隔的采样，控制器对采样信号进行数字处理并产生相应的离散控制信号，通过零阶保持器将控制信号输入轨道推进器使其产生连续的控制推力驱动航天器进行相应的轨道机动。这个过程实际上是一个采样控制过程，采样点的间隔时间是采样控制的重要参数，也可以把这个采样间隔时间看作是数字计算机的处理周期。

综上，目前采用连续信号形式设计控制器时通常假设测量信号和控制信号均为严格的实时信号，忽略了数字控制器的处理周期，因此在实际应用中难以获得预期的控制效果。此外，由于多种复杂因素的影响，航天器轨道推进器在采样时刻产生的推力与控制器计算的期望推力之间存在难以测定的偏差，这也将很大程度上影响轨道机动的精确性和安全性。

发明内容

本发明为解决采用现有的航天器相对运动的采样控制方法忽略了数字控制器的处理周期和偏差，影响航天器轨道的精确性和安全性的问题，进而提供了一种航天器相对运动的采样控制方法。

本发明为解决上述技术问题采取的技术方案是：所述航天器相对运动的采样控制方法由以下步骤实现的：

1.一种航天器相对运动的采样控制方法，其特征在于所述采样控制方法由以下步骤实现的：

步骤A、建立航天器相对运动动力学模型：

设两个航天器为追踪航天器和目标航天器，目标轨道为近似圆轨道，以目标航天器作为原点建立相对运动坐标系

将目标航天器的质心作为坐标系原点o，ox轴位于目标航天器轨道平面内，正向为地心指向航天器方向；oy轴为目标航天器运行方向；oz轴垂直于轨道平面并与其他两轴构成右手直角坐标系；

设追踪航天器相对于目标航天器的相对位置在x，y及z轴上的分量为x(t)、y(t)和z(t)，相对运动速度在相应坐标轴上的分量为

和

则相对运动状态向量为

设u_x(t)、u_y(t)和u_z(t)分别为作用在x、y和z轴上的控制推力，则控制输入向量定义为u(t)＝[u_x(t)，u_y(t)，u_z(t)]^T；追踪航天器质量为m，则相对运动的状态空间的系统方程可以写为：

$\stackrel{\·}{x}\left(t\right)=\mathrm{Ax}\left(t\right)+\mathrm{Bu}\left(t\right)---\left(1\right)$

式中A为系统状态矩阵，B为输入矩阵，分别有如下形式：

$A=\left[\begin{array}{cccccc}0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 1\\ 3n^2& 0& 0& 0& 2n& 0\\ 0& 0& 0& -2n& 0& 0\\ 0& 0& -n^2& 0& 0& 0\end{array}\right],$ $B=\frac{1}{m}\left[\begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ 0& 0& 0\\ 1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right]$

其中n为目标航天器的运行角速度；

步骤B、对两个航天器相对状态进行采样：

在追踪航天器和目标航天器相对运动过程中，采样器在采样时刻对追踪航天器和目标航天器相对状态进行采样，控制器根据采样信号计算此时刻的控制推力并产生离散形式的控制信号，控制信号通过零阶保持器驱动轨道推进器输出相应的连续控制推力；

步骤B1、设t_k是采样点时刻，两个航天器的相对运动状态在t₁、t₂、...、t_k、t_k+1、...时刻被控制器采集，处于t_k≤t＜t_k+1时间段内的状态均被认为是t_k时刻的状态进行处理；同样，与t_k时刻运动状态相对应的推力控制信号也以采样信号形式输出到零阶保持器，进而驱动推进器在t_k≤t＜t_k+1时段内对追踪航天器以此推力进行相应的机动控制；可见，对于t_k≤t＜t_k+1，系统方程(1)中连续形式的控制输入向量u(t)转化为采样点形式的控制输入向量u(t_k)，其形式为：

u(t_k)＝Kx(t_k)      (2)

两个航天器相对运动过程中有限推力条件由下式表示：

|u_i)t_k)|u_i，max，i＝x，y，z       (3)

其中u_i，max(i＝x，y，z)为x、y和z轴上的控制推力上界；

步骤B2、确定实际推力为u_r与控制律期望推力为u_d之间的关系

假设推进器产生的推力值同期望推力值之间的偏差分布在一个确定的范围内，设实际推力为u_r，控制律期望推力为u_d；

当控制律期望推力u_d＝0时，即推进器关闭，此时推进器的非线性特性也无从体现，因此输出值u_r＝0；但当控制律期望推力不为零时，推进器开始工作，随着推力需求增大，其非线性影响产生的推力偏差也相应增大，而且此偏差通常难以测得；

在期望推力值附近假想出的一个扇形区域，扇形区域的上下边界线分别为u_r＝σ_hu_d和u_r＝σ_lu_d，实际推力值分布于此扇形区域内，设推进器实际输出的推力向量为

S(u(t_k))＝[sec_x(u_x(t_k))，sec_y(u_y(t_k))，sec_z(u_z(t_k))]^T     (3)

式中sec_i(u_i(t_k))(i＝x，y，z)为x、y和z轴上的实际输出推力，满足如下关系

σ_liu_i(t_k)≤sec_i(u_i(t_k))≤σ_hiu_i(t_k)，i＝x，y，z       (4)

其中σ_li(i＝x，y，z)为x、y和z轴上推力扇形区域的下界比例系数，σ_hi(i＝x，y，z)为x、y和z轴上推力扇形区域的上界比例系数；由(1)、(2)、(3)式，相对运动系统方程可转化为如下形式：

$\stackrel{\·}{x}\left(t\right)=\mathrm{Ax}\left(t\right)+\mathrm{BS}\left(u\left(t_k\right)\right)---\left(5\right)$

步骤C、利用步骤B中所述的扇形区域的上下边界线构造M和N矩阵，

公式如下：

$M=\frac{1}{2}\mathrm{diag}\{({\mathrm{\σ}}_{\mathrm{lx}}+{\mathrm{\σ}}_{\mathrm{hx}}),({\mathrm{\σ}}_{\mathrm{ly}}+{\mathrm{\σ}}_{\mathrm{hy}}),({\mathrm{\σ}}_{\mathrm{lz}}+{\mathrm{\σ}}_{\mathrm{hz}})\}---\left(6\right)$

$N=\frac{1}{2}\mathrm{diag}\{({\mathrm{\σ}}_{\mathrm{hx}}-{\mathrm{\σ}}_{\mathrm{lx}}),({\mathrm{\σ}}_{\mathrm{hy}}-{\mathrm{\σ}}_{\mathrm{ly}}),({\mathrm{\σ}}_{\mathrm{hz}}-{\mathrm{\σ}}_{\mathrm{lz}})\}---\left(7\right)$

式中diag{}表示对角矩阵，定义向量

η(t_k)＝S(u(t_k))-Mu(t_k)        (8)

由公式(8)可得到实际输出控制推力S(u(t_k))，见公式(9)

S(u(t_k))＝η(t_k)+Mu(t_k)           (9)

步骤D、求得相应的状态反馈控制律：

设相邻两个采样点的时间间隔上界为h，即t_k+1-t_k≤h；定义d(t)＝t-t_k，则d(t)满足d(t)≤h，且采样点t_k可以写为t_k＝t-(t-t_k)＝t-d(t)，采样点时刻的状态向量可写为

x(t_k)＝x(t-d(t))     (10)

则由(2)和(10)可得相应的状态反馈控制律为

u(t_k)＝Kx(t_k)＝Kx(t-d(t))       (11)

将公式(9)和公式(11)代入公式(5)，可将相对运动系统方程进一步转化为：

$\stackrel{\·}{x}\left(t\right)=\mathrm{Ax}\left(t\right)+\mathrm{B\η}\left(t_k\right)+\mathrm{BMKx}(t-d\left(t\right))---\left(12\right)$

步骤E、引入两个正定对称矩阵P和Q并定义如下李亚普诺夫泛函

由状态向量x(t)的定义可知，x(t)由一个非零向量收敛到一个零向量即意味着两航天器相对位置和相对速度均为零，则系统方程(12)的渐进稳定也就意味着追踪航天器与目标航天器能够实现交会，为了保证相对运动系统方程(12)的渐进稳定性，引入两个正定对称矩阵P和Q并定义如下李亚普诺夫泛函

V(t)＝V₁(t)+V₂(t)        (13)

其中

V₁(t)＝x^T(t)Px(t)， $V_2\left(t\right)={\∫}_{-h}^0{\∫}_{t+\mathrm{\β}}^t{\stackrel{\·}{x}}^T\left(\mathrm{\α}\right)Q\stackrel{\·}{x}\left(\mathrm{\α}\right)\mathrm{d\αd\β}---\left(14\right)$

根据相应定义及矩阵不等式相关结论，可得

其中

矩阵Γ₁、Γ₂和Γ₃由以下各式给定

${\mathrm{\Γ}}_1=\left[\begin{array}{cc}\mathrm{\Π}& \mathrm{PBMK}\\ *& {\mathrm{\ϵ}}_1^{-1}{K}^T\mathrm{NNK}\end{array}\right],$ ${\mathrm{\Γ}}_3=\left[\begin{array}{c}I\\ -I\end{array}\right]Q\left[\begin{array}{cc}I& I\end{array}\right]$

${\mathrm{\Γ}}_2=\left[\begin{array}{c}\mathrm{QA}\\ \mathrm{QBMK}\end{array}\right]{\mathrm{\Θ}}^{-1}\left[\begin{array}{cc}\mathrm{QA}& \mathrm{QBMK}\end{array}\right]+{\mathrm{\ϵ}}_2\left[\begin{array}{c}0\\ \mathrm{NK}\end{array}\right]\left[\begin{array}{cc}0& \mathrm{NK}\end{array}\right],$ $\mathrm{\Θ}=Q-{\mathrm{\ϵ}}_2^{-1}{\mathrm{QBB}}^{T}Q$

根据公式(15)，如果矩阵K能够满足下式

Γ₁+hΓ₂-h^-1Γ₃＜0       (16)

那么

即系统方程(12)渐近稳定，从而航天器能够实现交会，因此，将(16)式作为控制律设计过程的一个约束条件；

步骤F、求得交会过程完成并且推力满足公式(3)上界约束条件：

步骤F1、将航天器轨道机动推进器的有限推力条件写为如下形式

$x^T(t-d\left(t\right))K^T{U}_i^T{U}_i\mathrm{Kx}(t-d\left(t\right))\≤u_{i,\max }^2---\left(17\right)$

则有限推力条件可由以下不等式条件满足

$\mathrm{\&rho;}{K}^T{U}_i^T{U}_{i}K<u_{i,\max }^{2}P---\left(18\right)$

其中ρ为一个给定常数满足V(0)＜ρ，其中V(0)为(13)式在初始条件下的取值；可见，利用(18)式结合上一步中得到的(16)式求得的控制矩阵K即可保证交会过程完成并且推力满足上界约束条件(3)；

步骤F2、对(16)和(18)式进行求解，通过矩阵不等式变换将两式进一步转化为以下两个矩阵不等式

$\left[\begin{array}{cc}-h^{-1}{\mathrm{\&epsiv;}}_{2}I& \widetilde{\mathrm{\&Phi;}}\\ *& \widetilde{\mathrm{\&Psi;}}\end{array}\right]<0---\left(19\right)$

$\left[\begin{array}{cc}-X& \sqrt{\mathrm{\&rho;}}{Y}^T{U}_i^T\\ *& -\mathrm{\&mu;I}\end{array}\right]<0---\left(20\right)$

式中X＝P^-1，Y＝KX，

μ为一给定正数并满足

相应矩阵具有如下形式

$\widetilde{\mathrm{\&Phi;}}=\left[\begin{array}{cccccc}{B}^T& 0& 0& 0& 0& 0\end{array}\right]$

$\widetilde{\mathrm{\&Psi;}}=\left[\begin{array}{cc}-h(\tilde{Q}-2X)& {\widetilde{\mathrm{\&Psi;}}}_{12}\\ *& {\widetilde{\mathrm{\&Psi;}}}_{22}\end{array}\right],$ ${\widetilde{\mathrm{\&Psi;}}}_{22}=\left[\begin{array}{cc}{\widetilde{\mathrm{\&gamma;}}}_1& {\widetilde{\mathrm{\&gamma;}}}_2\\ *& \widetilde{\mathrm{\&Lambda;}}+{\widetilde{\mathrm{\&gamma;}}}_3+{\widetilde{\mathrm{\&gamma;}}}_4\end{array}\right]$

${\widetilde{\mathrm{\&Psi;}}}_{12}=\left[\begin{array}{ccccc}\mathrm{hX}& \mathrm{hBMY}& 0& 0& 0\end{array}\right]$

${\widetilde{\mathrm{\&gamma;}}}_1=\mathrm{AX}+{\mathrm{XA}}^T-h^{-1}\tilde{Q},$ ${\widetilde{\mathrm{\&gamma;}}}_2=\left[\begin{array}{cccc}\mathrm{BMY}+h^{-1}\tilde{Q}& B& 0& 0\end{array}\right]$

${\widetilde{\mathrm{\&gamma;}}}_3={\left[\begin{array}{cccc}I& 0& 0& 0\end{array}\right]}^T\mathrm{NY}\left[\begin{array}{cccc}0& 0& I& 0\end{array}\right]$

${\widetilde{\mathrm{\&gamma;}}}_4={\left[\begin{array}{cccc}I& 0& 0& 0\end{array}\right]}^T\mathrm{NY}\left[\begin{array}{cccc}0& 0& 0& I\end{array}\right]$

$\widetilde{\mathrm{\&Lambda;}}=\mathrm{diag}\left\{\begin{array}{cccc}-h^{-1}\tilde{Q}& -{\mathrm{\&epsiv;}}_2^{-1}I,& -{\mathrm{\&epsiv;}}_{1}I,& -{\left({\mathrm{h\&epsiv;}}_2\right)}^{-1}I\end{array}\right\}$

如果给定u_max，则μ给定，(19)式和(20)式是关于X、Y和

的线性矩阵不等式；

步骤G、利用MATLAB软件中线性矩阵不等式(LMI)工具箱求可行解：

利用MATLAB软件中线性矩阵不等式(LMI)工具箱对于(19)和(20)式进行求解得到其可行解

利用算得的X和Y矩阵通过下式计算状态反馈增益矩阵K

K＝YX^-1          (21)

至此，即得到满足设计要求的航天器相对运动的状态反馈采样控制律

u(t_k)＝Kx(t_k)。

本发明具有以下有益效果：

本发明的采样控制方法考虑了数字控制器的处理周期，并且考虑了采样时刻产生的推力与控制器计算的期望推力之间存在难以测定的偏差，与现有的航天器相对运动的采样控制方法相比，本发明的采样控制方法的状态反馈控制律可以使两航天器在相应控制推力作用下实现交会，所需推力在运行过程中航天器通过实时状态在轨确定，并且所需推力均在允许推力范围内，大大提高了航天器轨道机动的精确性和安全性，本发明在实际应用中可以获得预期的控制效果。

附图说明

图1是本发明的航天器相对运动的采样控制方法的流程图，图2是本发明航天器相对运动坐标系示意图(其中O为地球质心，1为追踪航天器，2为目标航天器)，图3是航天器相对运动采样控制系统示意图，图4是实际推力u_r与期望推力u_d的关系图(其中

表示推力扇形区域上界，

表示期望推力，

表示推力扇形区域下界)，图5是航天器相对位置在x轴和y轴上分量随时间变化曲线，图6是交会过程控制律期望推力随时间变化图，图7是交会过程实际控制推力随时间变化图，图8是δ＝0.05时对应的轨道机动过程所需的最大推力随采样时间变化图，图9是δ＝0.1时对应的轨道机动过程所需的最大推力随采样时间变化图，图10是δ＝0.15时对应的轨道机动过程所需的最大推力随采样时间变化图，图11是δ＝0.2时对应的轨道机动过程所需的最大推力随采样时间变化图。

具体实施方式

具体实施方式一：本实施方式的航天器相对运动的采样控制方法是由以下步骤实现的：

步骤A、建立航天器相对运动动力学模型：

设两个航天器为追踪航天器和目标航天器，目标轨道为近似圆轨道，以目标航天器作为原点建立相对运动坐标系(如图2所示)

将目标航天器的质心作为坐标系原点o，ox轴位于目标航天器轨道平面内，正向为地心指向航天器方向；oy轴为目标航天器运行方向；oz轴垂直于轨道平面并与其他两轴构成右手直角坐标系；

设追踪航天器相对于目标航天器的相对位置在x，y及z轴上的分量为x(t)、y(t)和z(t)，相对运动速度在相应坐标轴上的分量为

和

则相对运动状态向量为

设u_x(t)、u_y(t)和u_z(t)分别为作用在x、y和z轴上的控制推力，则控制输入向量定义为u(t)＝[u_x(t)，u_y(t)，y_z(t)]^T；追踪航天器质量为m，则相对运动的状态空间的系统方程可以写为：

$\stackrel{\&CenterDot;}{x}\left(t\right)=\mathrm{Ax}\left(t\right)+\mathrm{Bu}\left(t\right)---\left(1\right)$

式中A为系统状态矩阵，B为输入矩阵，分别有如下形式：

$A=\left[\begin{array}{cccccc}0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 1\\ 3n^2& 0& 0& 0& 2n& 0\\ 0& 0& 0& -2n& 0& 0\\ 0& 0& -n^2& 0& 0& 0\end{array}\right],$ $B=\frac{1}{m}\left[\begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ 0& 0& 0\\ 1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right]$

其中n为目标航天器的运行角速度；

步骤B、对两个航天器相对状态进行采样：

在追踪航天器和目标航天器相对运动过程中，采样器在采样时刻对追踪航天器和目标航天器相对状态进行采样，控制器根据采样信号计算此时刻的控制推力并产生离散形式的控制信号，控制信号通过零阶保持器驱动轨道推进器输出相应的连续控制推力(如图3所示)

步骤B1、设t_k是采样点时刻，两个航天器的相对运动状态在t₁、t₂、...、t_k、t_k+1、...时刻被控制器采集，处于t_k≤t＜t_k+1时间段内的状态均被认为是t_k时刻的状态进行处理；同样，与t_k时刻运动状态相对应的推力控制信号也以采样信号形式输出到零阶保持器，进而驱动推进器在t_k≤t＜t_k+1时段内对追踪航天器以此推力进行相应的机动控制；可见，对于t_k≤t＜t_k+1，系统方程(1)中连续形式的控制输入向量u(t)转化为采样点形式的控制输入向量u(t_k)，其形式为：

u(t_k)＝Kx(t_k)       (2)

两个航天器相对运动过程中有限推力条件由下式表示：

|u_i(t_k)|u_i，max，i＝x，y，z      (3)

其中uX(i＝x，y，z)为x、y和z轴上的控制推力上界；

步骤B2、确定实际推力为u_r，控制律期望推力为u_d之间的关系

在每个采样点时刻，由于推进器内部燃料损耗、摩擦等因素，推进器产生的推力值同期望推力值必然存在难以测定的偏差；假设推进器产生的推力值同期望推力值之间的偏差分布在一个确定的范围内，设实际推力为u_r，控制律期望推力为u_d，两个推力值之间的关系如图4所示；

当控制律期望推力u_d＝0时，即推进器关闭，此时推进器的非线性特性也无从体现，因此输出值u_r＝0；但当控制律期望推力不为零时，推进器开始工作，随着推力需求增大，其非线性影响产生的推力偏差也相应增大，而且此偏差通常难以测得；

在期望推力值附近假想出的一个扇形区域(如图4所示)，扇形区域的上下边界线分别为u_r＝σ_hu_d和u_r＝σ_lu_d，实际推力值分布于此扇形区域内，设推进器实际输出的推力向量为

S(u(t_k))＝[sec_x(u_x(t_k))，sec_y(u_y(t_k))，sec_z(u_z(t_k))]^T     (3)

式中sec_i(u_i(t_k))(i＝x，y，z)为x、y和z轴上的实际输出推力，满足如下关系

σ_liu_i(t_k)≤sec_i(u_i(t_k))≤σ_hiu_i(t_k)，i＝x，y，z      (4)

其中σ_li(i＝x，y，z)为x、y和z轴上推力扇形区域的下界比例系数，σ_hi(i＝x，y，z)为x、y和z轴上推力扇形区域的上界比例系数；由(1)、(2)、(3)式，相对运动系统方程可转化为如下形式：

$\stackrel{\&CenterDot;}{x}\left(t\right)=\mathrm{Ax}\left(t\right)+\mathrm{BS}\left(u\left(t_k\right)\right)---\left(5\right)$

步骤C、利用步骤B中所述的扇形区域的上下边界线构造M和N矩阵，公式如下：

$M=\frac{1}{2}\mathrm{diag}\{({\mathrm{\&sigma;}}_{\mathrm{lx}}+{\mathrm{\&sigma;}}_{\mathrm{hx}}),({\mathrm{\&sigma;}}_{\mathrm{ly}}+{\mathrm{\&sigma;}}_{\mathrm{hy}}),({\mathrm{\&sigma;}}_{\mathrm{lz}}+{\mathrm{\&sigma;}}_{\mathrm{hz}})\}---\left(6\right)$

$N=\frac{1}{2}\mathrm{diag}\{({\mathrm{\&sigma;}}_{\mathrm{hx}}-{\mathrm{\&sigma;}}_{\mathrm{lx}}),({\mathrm{\&sigma;}}_{\mathrm{hy}}-{\mathrm{\&sigma;}}_{\mathrm{ly}}),({\mathrm{\&sigma;}}_{\mathrm{hz}}-{\mathrm{\&sigma;}}_{\mathrm{lz}})\}---\left(7\right)$

式中diag{}表示对角矩阵，定义向量

η(t_k)＝S(u(t_k))-Mu(t_k)        (8)

由公式(8)可得到实际输出控制推力S(u(t_k))，见公式(9)

S(u(t_k))＝η(t_k)+Mu(t_k)         (9)

步骤D、求得相应的状态反馈控制律：

设相邻两个采样点的时间间隔上界为h，即t_k+1-t_k≤h；定义d(t)＝t-t_k，则d(t)满足d(t)≤h，且采样点t_k可以写为t_k＝t-(t-t_k)＝t-d(t)，采样点时刻的状态向量可写为

x(t_k)＝x(t-d(t))      (10)

则由(2)和(10)可得相应的状态反馈控制律为

u(t_k)＝Kx(t_k)＝Kx(t-d(t))        (11)

将公式(9)和公式(11)代入公式(5)，可将相对运动系统方程进一步转化为：

$\stackrel{\&CenterDot;}{x}\left(t\right)=\mathrm{Ax}\left(t\right)+\mathrm{B\&eta;}\left(t_k\right)+\mathrm{BMKx}(t-d\left(t\right))---\left(12\right)$

步骤E、引入两个正定对称矩阵P和Q并定义如下李亚普诺夫泛函

由状态向量x(t)的定义可知，x(t)由一个非零向量收敛到一个零向量即意味着两航天器相对位置和相对速度均为零，则系统方程(12)的渐进稳定也就意味着追踪航天器与目标航天器能够实现交会，为了保证相对运动系统方程(12)的渐进稳定性，引入两个正定对称矩阵P和Q并定义如下李亚普诺夫泛函

V(t)＝V₁(t)+V₂(t)     (13)

其中

V₁(t)＝x^T(t)Px(t)， $V_2\left(t\right)={\&Integral;}_{-h}^0{\&Integral;}_{t+\mathrm{\&beta;}}^t{\stackrel{\&CenterDot;}{x}}^T\left(\mathrm{\&alpha;}\right)Q\stackrel{\&CenterDot;}{x}\left(\mathrm{\&alpha;}\right)\mathrm{d\&alpha;d\&beta;}---\left(14\right)$

根据相应定义及矩阵不等式相关结论，可得

其中

矩阵Γ₁、Γ₂和Γ₃由以下各式给定

${\mathrm{\&Gamma;}}_1=\left[\begin{array}{cc}\mathrm{\&Pi;}& \mathrm{PBMK}\\ *& {\mathrm{\&epsiv;}}_1^{-1}{K}^T\mathrm{NNK}\end{array}\right],$ ${\mathrm{\&Gamma;}}_3=\left[\begin{array}{c}I\\ -I\end{array}\right]Q\left[\begin{array}{cc}I& I\end{array}\right]$

${\mathrm{\&Gamma;}}_2=\left[\begin{array}{c}\mathrm{QA}\\ \mathrm{QBMK}\end{array}\right]{\mathrm{\&Theta;}}^{-1}\left[\begin{array}{cc}\mathrm{QA}& \mathrm{QBMK}\end{array}\right]+{\mathrm{\&epsiv;}}_2\left[\begin{array}{c}0\\ \mathrm{NK}\end{array}\right]\left[\begin{array}{cc}0& \mathrm{NK}\end{array}\right],$ $\mathrm{\&Theta;}=Q-{\mathrm{\&epsiv;}}_2^{-1}{\mathrm{QBB}}^{T}Q$

根据公式(15)，如果矩阵K能够满足下式

Γ₁+hΓ₂-h^-1Γ₃＜0        (16)

那么即系统方程(12)渐近稳定，从而航天器能够实现交会，因此，将(16)式作为控制律设计过程的一个约束条件；

步骤F、保证交会过程完成并且推力满足公式(3)上界约束条件：

步骤F1、将航天器轨道机动推进器的有限推力条件写为如下形式

$x^T(t-d\left(t\right))K^T{U}_i^T{U}_i\mathrm{Kx}(t-d\left(t\right))\&le;u_{i,\max }^2---\left(17\right)$

则有限推力条件可由以下不等式条件满足

$\mathrm{\&rho;}{K}^T{U}_i^T{U}_{i}K<u_{i,\max }^{2}P---\left(18\right)$

其中ρ为一个给定常数满足V(0)＜ρ，其中V(0)为(13)式在初始条件下的取值；可见，利用(18)式结合上一步中得到的(16)式求得的控制矩阵K即可保证交会过程完成并且推力满足上界约束条件(3)；

步骤F2、对(16)和(18)式进行求解，通过矩阵不等式变换将两式进一步转化为以下两个矩阵不等式

$\left[\begin{array}{cc}-h^{-1}{\mathrm{\&epsiv;}}_{2}I& \widetilde{\mathrm{\&Phi;}}\\ *& \widetilde{\mathrm{\&Psi;}}\end{array}\right]<0---\left(19\right)$

$\left[\begin{array}{cc}-X& \sqrt{\mathrm{\&rho;}}{Y}^T{U}_i^T\\ *& -\mathrm{\&mu;I}\end{array}\right]<0---\left(20\right)$

式中X＝P^-1，Y＝KX，μ为一给定正数并满足

相应矩阵具有如下形式

$\widetilde{\mathrm{\&Phi;}}=\left[\begin{array}{cccccc}{B}^T& 0& 0& 0& 0& 0\end{array}\right]$

$\widetilde{\mathrm{\&Psi;}}=\left[\begin{array}{cc}-h(\tilde{Q}-2X)& {\widetilde{\mathrm{\&Psi;}}}_{12}\\ *& {\widetilde{\mathrm{\&Psi;}}}_{22}\end{array}\right],$ ${\widetilde{\mathrm{\&Psi;}}}_{22}=\left[\begin{array}{cc}{\widetilde{\mathrm{\&gamma;}}}_1& {\widetilde{\mathrm{\&gamma;}}}_2\\ *& \widetilde{\mathrm{\&Lambda;}}+{\widetilde{\mathrm{\&gamma;}}}_3+{\widetilde{\mathrm{\&gamma;}}}_4\end{array}\right]$

${\widetilde{\mathrm{\&Psi;}}}_{12}=\left[\begin{array}{ccccc}\mathrm{hX}& \mathrm{hBMY}& 0& 0& 0\end{array}\right]$

${\widetilde{\mathrm{\&gamma;}}}_1=\mathrm{AX}+{\mathrm{XA}}^T-h^{-1}\tilde{Q},$ ${\widetilde{\mathrm{\&gamma;}}}_2=\left[\begin{array}{cccc}\mathrm{BMY}+h^{-1}\tilde{Q}& B& 0& 0\end{array}\right]$

${\widetilde{\mathrm{\&gamma;}}}_3={\left[\begin{array}{cccc}I& 0& 0& 0\end{array}\right]}^T\mathrm{NY}\left[\begin{array}{cccc}0& 0& I& 0\end{array}\right]$

${\widetilde{\mathrm{\&gamma;}}}_4={\left[\begin{array}{cccc}I& 0& 0& 0\end{array}\right]}^T\mathrm{NY}\left[\begin{array}{cccc}0& 0& 0& I\end{array}\right]$

$\widetilde{\mathrm{\&Lambda;}}=\mathrm{diag}\left\{\begin{array}{cccc}-h^{-1}\tilde{Q}& -{\mathrm{\&epsiv;}}_2^{-1}I,& -{\mathrm{\&epsiv;}}_{1}I,& -{\left({\mathrm{h\&epsiv;}}_2\right)}^{-1}I\end{array}\right\}$

如果给定u_max，则μ给定，(19)式和(20)式是关于X、Y和

的线性矩阵不等式；

步骤G、利用MATLAB软件中线性矩阵不等式(LMI)工具箱求可行解：

利用MATLAB软件中线性矩阵不等式(LMI)工具箱对于(19)和(20)式进行求解得到其可行解

利用算得的X和Y矩阵通过下式计算状态反馈增益矩阵K

K＝YX^-1       (21)

至此，即得到满足设计要求的航天器相对运动的状态反馈采样控制律

u(t_k)＝Kx(t_k)。

本发明方法的实例验证：

1)目标航天器质量：200kg

2)目标航天器运行轨道半径：42241km

3)目标航天器轨道运行平均角速度：0.001117rad/s

4)初始时刻两航天器的相对状态：[100，150，0，0，0，0]

5)设定脉冲推力上界为500N

6)追踪航天器轨道推进器产生的实际推力与期望推力满足如下关系

S(u(t_k))＝u(t_k)+δu(t_k)sin[u(t_k)]

式中δ为一任意给定的常数值，则非线性影响程度可通过调整参数δ的大小来实现；

控制律求解：

设定采样间隔时间上界h＝0.1s，推力非线性程度δ＝0.1，依据上述计算过程，利用MATLAB软件线性矩阵不等式(LMI)工具箱对不等式(19)式和(20)式进行求解，得到状态反馈增益矩阵K为如下形式

$K=\left[\begin{array}{cccccc}-20.9871& -1.3618& -4.0628& -84.2685& 7.3177& -0.5152\\ -2.8031& -17.9294& -3.5836& -3.1000& -63.1909& 0.6210\\ -3.0069& -0.2291& -22.6558& -3.2550& 18.9533& -84.6244\end{array}\right]$

控制律作用效果：

根据上述结果，得到状态反馈脉冲控制律u(t_k)＝Kx(t_k)，将此控制律应用于追踪航天器，使其从初始位置开始自主确定交会过程所需控制推力进行运动(如图5所示)

采用设计的控制律追踪航天器可以在轨根据实时相对运动状态自主计算所需控制推力的大小，以y-轴为例，如图6和图7所示的交会过程中y-轴推力变化情况以及实际推力和期望推力间的对比图，由图6和图7可见，控制推力符合给定的有限脉冲推力条件，并且在推力存在偏差情况下航天器依然能够实现交会；

此外，对于反馈增益矩阵K的求解，参数h和δ的取值具有一定的限度：当采样间隔过大或推进器非线性影响较大时，可能会造成线性矩阵不等式无法求解；因此，此控制律设计方法对于h和δ的容忍程度也是对其进行评价的重要指标；表1列出对于不同的h取值，能够保证控制器K存在的最大δ值，表2列出不同δ情况下，保证控制器K存在的最大采样间隔上界h_max，表1和表2如下：

表1 不同采样点间隔上界h对应的最大可容忍非线性程度δ

表2 不同非线性程度δ对应的最大可容忍采样间隔上界h

参数h和δ的不同取值在求解控制律或对轨道机动过程进行分析时都具有非常重要的影响，图8～11给出了δ＝0.05、δ＝0.1、δ＝0.15和δ＝0.2四种情况下不同h值对应的轨道机动过程所需的最大推力；

综合以上各图可见，对于允许的采样间隔时间上界h和推力非线性程度参数δ，应用所设计的状态反馈控制律可以使两航天器在相应控制推力作用下实现交会，所需推力在运行过程中航天器通过实时状态在轨确定，并且所需推力均在允许推力范围内。

Claims (1)

1.一种航天器相对运动的采样控制方法，其特征在于所述采样控制方法由以下步骤实现的：

步骤A、建立航天器相对运动动力学模型：

设两个航天器为追踪航天器和目标航天器，目标轨道为近似圆轨道，以目标航天器作为原点建立相对运动坐标系

将目标航天器的质心作为坐标系原点o，ox轴位于目标航天器轨道平面内，正向为地心指向航天器方向；oy轴为目标航天器运行方向；oz轴垂直于轨道平面并与其他两轴构成右手直角坐标系；

设追踪航天器相对于目标航天器的相对位置在x，y及z轴上的分量为x(t)、y(t)和z(t)，相对运动速度在相应坐标轴上的分量为

和则相对运动状态向量为

设u_x(t)、u_y(t)和u_z(t)分别为作用在x、y和z轴上的控制推力，则控制输入向量定义为u(t)＝[u_x(t)，u_y(t)，u_z(t)]^T；追踪航天器质量为m，则相对运动的状态空间的系统方程可以写为：

$\stackrel{\&CenterDot;}{x}\left(t\right)=\mathrm{Ax}\left(t\right)+\mathrm{Bu}\left(t\right)---\left(1\right)$

式中A为系统状态矩阵，B为输入矩阵，分别有如下形式：

$A=\left[\begin{array}{cccccc}0& 0& 0& 1& 0& 0\\ 0& 0& 0& 0& 1& 0\\ 0& 0& 0& 0& 0& 1\\ 3n^2& 0& 0& 0& 2n& 0\\ 0& 0& 0& -2n& 0& 0\\ 0& 0& -n^2& 0& 0& 0\end{array}\right],$ $B=\frac{1}{m}\left[\begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ 0& 0& 0\\ 1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right]$

其中n为目标航天器的运行角速度；

步骤B、对两个航天器相对状态进行采样：

在追踪航天器和目标航天器相对运动过程中，采样器在采样时刻对追踪航天器和目标航天器相对状态进行采样，控制器根据采样信号计算此时刻的控制推力并产生离散形式的控制信号，控制信号通过零阶保持器驱动轨道推进器输出相应的连续控制推力；

步骤B1、设t_k是采样点时刻，两个航天器的相对运动状态在t₁、t₂、...、t_k、t_k+1、...时刻被控制器采集，处于t_k≤t＜t_k+1时间段内的状态均被认为是t_k时刻的状态进行处理；同样，与t_k时刻运动状态相对应的推力控制信号也以采样信号形式输出到零阶保持器，进而驱动推进器在t_k≤t＜t_k+1时段内对追踪航天器以此推力进行相应的机动控制；可见，对于t_k≤t＜t_k+1，系统方程(1)中连续形式的控制输入向量u(t)转化为采样点形式的控制输入向量u(t_k)，其形式为：

u(t_k)＝Kx(t_k)       (2)

两个航天器相对运动过程中有限推力条件由下式表示：

|u_i(t_k)|≤u_i，max，i＝x，y，z      (3)

其中u_i，max(i＝x，y，z)为x、y和z轴上的控制推力上界；

步骤B2、确定实际推力为u_r与控制律期望推力为u_d之间的关系

假设推进器产生的推力值同期望推力值之间的偏差分布在一个确定的范围内，设实际推力为u_r，控制律期望推力为u_d；

当控制律期望推力u_d＝0时，即推进器关闭，此时推进器的非线性特性也无从体现，因此输出值u_r＝0；但当控制律期望推力不为零时，推进器开始工作，随着推力需求增大，其非线性影响产生的推力偏差也相应增大，而且此偏差通常难以测得；

在期望推力值附近假想出的一个扇形区域，扇形区域的上下边界线分别为u_r＝σ_hu_d和u_r＝σ_lu_d，实际推力值分布于此扇形区域内，设推进器实际输出的推力向量为

S(u(t_k))＝[sec_x(u_x(t_k))，sec_y(u_y(t_k))，sec_z(u_z(t_k))]^T      (3)

式中sec_i(u_i(t_k))(i＝x，y，z)为x、y和z轴上的实际输出推力，满足如下关系

σ_liu_i(t_k)≤sec_i(u_i(t_k))≤σ_hiu_i(t_k)，i＝x，y，z       (4)

其中σ_li(i＝x，y，z)为x、y和z轴上推力扇形区域的下界比例系数，σ_hi(i＝x，y，z)为x、y和z轴上推力扇形区域的上界比例系数；由(1)、(2)、(3)式，相对运动系统方程可转化为如下形式：

$\stackrel{\&CenterDot;}{x}\left(t\right)=\mathrm{Ax}\left(t\right)+\mathrm{BS}\left(u\left(t_k\right)\right)---\left(5\right)$

步骤C、利用步骤B中所述的扇形区域的上下边界线构造M和N矩阵，公式如下：

$M=\frac{1}{2}\mathrm{diag}\{({\mathrm{\&sigma;}}_{\mathrm{lx}}+{\mathrm{\&sigma;}}_{\mathrm{hx}}),({\mathrm{\&sigma;}}_{\mathrm{ly}}+{\mathrm{\&sigma;}}_{\mathrm{hy}}),({\mathrm{\&sigma;}}_{\mathrm{lz}}+{\mathrm{\&sigma;}}_{\mathrm{hz}})\}---\left(6\right)$

$N=\frac{1}{2}\mathrm{diag}\{({\mathrm{\&sigma;}}_{\mathrm{hx}}-{\mathrm{\&sigma;}}_{\mathrm{lx}}),({\mathrm{\&sigma;}}_{\mathrm{hy}}-{\mathrm{\&sigma;}}_{\mathrm{ly}}),({\mathrm{\&sigma;}}_{\mathrm{hz}}-{\mathrm{\&sigma;}}_{\mathrm{lz}})\}---\left(7\right)$

式中diag{}表示对角矩阵，定义向量

η(t_k)＝S(u(t_k))-Mu(t_k)        (8)

由公式(8)可得到实际输出控制推力S(u(t_k))，见公式(9)

S(u(t_k))＝η(t_k)+Mu(t_k)       (9)

步骤D、求得相应的状态反馈控制律：

设相邻两个采样点的时间间隔上界为h，即t_k+1-t_k≤h；定义d(t)＝t-t_k，则d(t)满足d(t)≤h，且采样点t_k可以写为t_k＝t-(t-t_k)＝t-d(t)，采样点时刻的状态向量可写为

x(t_k)＝x(t-d(t))        (10)

则由(2)和(10)可得相应的状态反馈控制律为

u(t_k)＝Kx(t_k)＝Kx(t-d(t))      (11)

将公式(9)和公式(11)代入公式(5)，可将相对运动系统方程进一步转化为：

$\stackrel{\&CenterDot;}{x}\left(t\right)=\mathrm{Ax}\left(t\right)+\mathrm{B\&eta;}\left(t_k\right)+\mathrm{BMKx}(t-d\left(t\right))---\left(12\right)$

步骤E、引入两个正定对称矩阵P和Q并定义如下李亚普诺夫泛函

由状态向量x(t)的定义可知，x(t)由一个非零向量收敛到一个零向量即意味着两航天器相对位置和相对速度均为零，则系统方程(12)的渐进稳定也就意味着追踪航天器与目标航天器能够实现交会，为了保证相对运动系统方程(12)的渐进稳定性，引入两个正定对称矩阵P和Q并定义如下李亚普诺夫泛函

V(t)＝V₁(t)+V₂(t)        (13)

其中

V₁(t)＝x^T(t)Px(t)， $V_2\left(t\right)={\&Integral;}_{-h}^0{\&Integral;}_{t+\mathrm{\&beta;}}^t{\stackrel{\&CenterDot;}{x}}^T\left(\mathrm{\&alpha;}\right)Q\stackrel{\&CenterDot;}{x}\left(\mathrm{\&alpha;}\right)\mathrm{d\&alpha;d\&beta;}---\left(14\right)$

根据相应定义及矩阵不等式相关结论，可得

其中矩阵Γ₁、Γ₂和Γ₃由以下各式给定

${\mathrm{\&Gamma;}}_1=\left[\begin{array}{cc}\mathrm{\&Pi;}& \mathrm{PBMK}\\ *& {\mathrm{\&epsiv;}}_1^{-1}{K}^T\mathrm{NNK}\end{array}\right],$ ${\mathrm{\&Gamma;}}_3=\left[\begin{array}{c}I\\ -I\end{array}\right]Q\left[\begin{array}{cc}I& I\end{array}\right]$

${\mathrm{\&Gamma;}}_2=\left[\begin{array}{c}\mathrm{QA}\\ \mathrm{QBMK}\end{array}\right]{\mathrm{\&Theta;}}^{-1}\left[\begin{array}{cc}\mathrm{QA}& \mathrm{QBMK}\end{array}\right]+{\mathrm{\&epsiv;}}_2\left[\begin{array}{c}0\\ \mathrm{NK}\end{array}\right]\left[\begin{array}{cc}0& \mathrm{NK}\end{array}\right],$ $\mathrm{\&Theta;}=Q-{\mathrm{\&epsiv;}}_2^{-1}{\mathrm{QBB}}^{T}Q$

根据公式(15)，如果矩阵K能够满足下式

Γ₁+hΓ₂-h^-1Γ₃＜0      (16)

那么即系统方程(12)渐近稳定，从而航天器能够实现交会，因此，将(16)式作为控制律设计过程的一个约束条件；

步骤F、求得交会过程完成并且推力满足公式(3)上界约束条件：

步骤F1、将航天器轨道机动推进器的有限推力条件写为如下形式

$x^T(t-d\left(t\right))K^T{U}_i^T{U}_i\mathrm{Kx}(t-d\left(t\right))\&le;u_{i,\max }^2---\left(17\right)$

则有限推力条件可由以下不等式条件满足

$\mathrm{\&rho;}{K}^T{U}_i^T{U}_{i}K<u_{i,\max }^{2}P---\left(18\right)$

其中ρ为一个给定常数满足V(0)＜ρ，其中V(0)为(13)式在初始条件下的取值；可见，利用(18)式结合上一步中得到的(16)式求得的控制矩阵K即可保证交会过程完成并且推力满足上界约束条件(3)；

步骤F 2、对(16)和(18)式进行求解，通过矩阵不等式变换将两式进一步转化为以下两个矩阵不等式

$\left[\begin{array}{cc}-h^{-1}{\mathrm{\&epsiv;}}_{2}I& \widetilde{\mathrm{\&Phi;}}\\ *& \widetilde{\mathrm{\&Psi;}}\end{array}\right]<0---\left(19\right)$

$\left[\begin{array}{cc}-X& \sqrt{\mathrm{\&rho;}}{Y}^T{U}_i^T\\ *& -\mathrm{\&mu;I}\end{array}\right]<0---\left(20\right)$

式中X＝P^-1，Y＝KX，

μ为一给定正数并满足

相应矩阵具有如下形式

$\widetilde{\mathrm{\&Phi;}}=\left[\begin{array}{cccccc}{B}^T& 0& 0& 0& 0& 0\end{array}\right]$

$\widetilde{\mathrm{\&Psi;}}=\left[\begin{array}{cc}-h(\tilde{Q}-2X)& {\widetilde{\mathrm{\&Psi;}}}_{12}\\ *& {\widetilde{\mathrm{\&Psi;}}}_{22}\end{array}\right],$ ${\widetilde{\mathrm{\&Psi;}}}_{22}=\left[\begin{array}{cc}{\widetilde{\mathrm{\&gamma;}}}_1& {\widetilde{\mathrm{\&gamma;}}}_2\\ *& \widetilde{\mathrm{\&Lambda;}}+{\widetilde{\mathrm{\&gamma;}}}_3+{\widetilde{\mathrm{\&gamma;}}}_4\end{array}\right]$

${\widetilde{\mathrm{\&Psi;}}}_{12}=\left[\begin{array}{ccccc}\mathrm{hX}& \mathrm{hBMY}& 0& 0& 0\end{array}\right]$

${\widetilde{\mathrm{\&gamma;}}}_1=\mathrm{AX}+{\mathrm{XA}}^T-h^{-1}\tilde{Q},$ ${\widetilde{\mathrm{\&gamma;}}}_2=\left[\begin{array}{cccc}\mathrm{BMY}+h^{-1}\tilde{Q}& B& 0& 0\end{array}\right]$

${\widetilde{\mathrm{\&gamma;}}}_3={\left[\begin{array}{cccc}I& 0& 0& 0\end{array}\right]}^T\mathrm{NY}\left[\begin{array}{cccc}0& 0& I& 0\end{array}\right]$

${\widetilde{\mathrm{\&gamma;}}}_4={\left[\begin{array}{cccc}I& 0& 0& 0\end{array}\right]}^T\mathrm{NY}\left[\begin{array}{cccc}0& 0& 0& I\end{array}\right]$

$\widetilde{\mathrm{\&Lambda;}}=\mathrm{diag}\left\{\begin{array}{cccc}-h^{-1}\tilde{Q}& -{\mathrm{\&epsiv;}}_2^{-1}I,& -{\mathrm{\&epsiv;}}_{1}I,& -{\left({\mathrm{h\&epsiv;}}_2\right)}^{-1}I\end{array}\right\}$

如果给定u_max，则μ给定，(19)式和(20)式是关于X、Y和

的线性矩阵不等式；

步骤G、利用MATLAB软件中线性矩阵不等式(LMI)工具箱求可行解：

利用MATLAB软件中线性矩阵不等式(LMI)工具箱对于(19)和(20)式进行求解得到其可行解

利用算得的X和Y矩阵通过下式计算状态反馈增益矩阵K

K＝YX^-1          (21)

至此，即得到满足设计要求的航天器相对运动的状态反馈采样控制律

u(t_k)＝Kx(t_k)。

Images