CN106379555A - 一种考虑j2摄动的近地卫星有限推力最优变轨方法 - Google Patents

Abstract

本发明公开的一种考虑J2摄动的近地卫星有限推力最优变轨方法，涉及近地卫星最优变轨方法，属于航空航天领域。本发明建立考虑J2摄动的航天器动力学方程；根据始末航天器状态，求解燃耗最优性能指标所对应的哈密顿函数H；根据庞特里亚金极大值原理求解最优控制率以及相应协态方程；根据变分法最优控制条件构造满足燃耗最优的两点边值问题打靶方程组；根据航天器始末状态和飞行时间，对两点边值问题打靶方程组进行求解，得到由最优推力控制方向和发动机开关函数构成的最优控制率；利用所述的控制率对卫星变轨进行控制，实现提高近地卫星转移过程的实时性和变轨精度，降低变轨转移所需的燃料消耗。本发明适用于发动机为有限推力模型下近地卫星变轨过程。

Description

一种考虑J2摄动的近地卫星有限推力最优变轨方法

技术领域

本发明是一种考虑J2摄动的近地卫星最优变轨方法，适用于发动机为有限推力模型下近地卫星变轨过程，属于航空航天技术领域。

背景技术

卫星变轨是航天器在太空运行工作中最常见的操作之一。相比于火箭发动机的大推力矢量，航天器上的发动机一般推力较小，在实施推力的过程中，单次推力时间和推力弧段相对较长，导致推力矢量在惯性空间是一个实时变化量，仍然采用脉冲设计并不能合理的反映这种控制过程。因此，对于有限推力的变轨策略要采用实时控制设计。而对于在轨航天器，其燃料质量决定了航天器的工作寿命，因此在变轨控制设计过程中，考虑燃耗最优的变轨设计策略能极大的降低变轨成本，从而增加航天器的在轨寿命。同时在考虑近地变轨过程中，J2摄动成为影响变轨结果的最主要因素，引入J2摄动影响可以较好的保证变轨精度。

在已发展的关于有限推力变轨方法中在先技术[1](参见Haberkorn T,MartinonP,Gergaud J.Low Thrust Minimum-Fuel Orbital Transfer:A Homotopic Approach[J].Journal of Guidance Control&Dynamics,2004,27(6):1046-1060.)给出了间接法的有限推力变轨优化策略，该方法采用同伦逼近的策略，通过变分法和庞特里亚金极大值原理推导出了满足最优控制策略的变轨控制率。该方法从理论上严格满足变轨的最优性条件，但在实际求解过程中为满足其一阶必要条件需要对其协态变量的初值进行较为正确的估计，但协态变量的物理意义远不如状态变量的物理意义明确，因此该方法对于初值猜测问题具有较高的要求。在求解过程中由于初值猜测的限制无法得到满足条件的可行解。

在先技术[2](参见Jiang F,Baoyin H,Li J.Practical Techniques for Low-Thrust Trajectory Optimization with Homotopic Approach[J].Journal of GuidanceControl&Dynamics,2012,1(35):245-258.)给出有限推力变轨方法的进一步改进，通过对协态的归一化处理，即通过对原性能指标乘以任意正数得到增广协态。该方法并不影响结果的最优性，同时极大地缩小了协态变量的可搜索范围，从而提高了初始猜测的优化效率。但该技术的研究仍简单的停留在二体引力模型上，只适用于不考虑摄动的变轨转移问题，对于精度要求较高的近地转移变轨，需要考虑J2摄动等因素时则无法解决。

发明内容

本发明公开的一种考虑J2摄动的近地卫星有限推力最优变轨方法要解决的技术问题是，提高近地卫星转移过程的实时性和变轨精度，降低变轨转移所需的燃料消耗。

本发明的目的是通过下述技术方案实现的：

本发明为公开的一种考虑J2摄动的近地卫星有限推力最优变轨方法，通过在地心赤道惯性坐标系中建立考虑J2摄动的航天器动力学方程。根据始末航天器状态，求解燃耗最优性能指标所对应的哈密顿函数H。根据庞特里亚金极大值原理求解最优控制率以及相应的协态方程。根据变分法最优控制条件构造满足燃耗最优的两点边值问题打靶方程组。根据航天器始末状态和飞行时间，对两点边值问题打靶方程组进行求解，得到由最优推力控制方向和发动机开关函数构成的最优控制率。利用所述的控制率对卫星变轨进行控制，实现提高近地卫星转移过程的实时性和变轨精度，降低变轨转移所需的燃料消耗。

本发明为公开的一种考虑J2摄动的近地卫星有限推力最优变轨方法，包括如下步骤：

步骤一：在地心赤道惯性坐标系中建立考虑J2摄动的航天器动力学方程。所述的J2摄动为非球形引力摄动。

首先考虑引入发动机的开关控制量u，在无摄动条件下笛卡尔坐标系中的有限推力航天器动力学方程公式(1)所示，此时笛卡尔坐标系选为地心赤道惯性坐标系，x轴指向春分点方向，z轴与地球自转轴一致，且指向北，y轴与x轴、z轴满足右手定则。

$\begin{array}{c}\stackrel{\·}{r}=v\\ \stackrel{\·}{v}=-\frac{\mathrm{\μ}}{r^3}r+\frac{T_{max}u}{m}\mathrm{\α}\\ \stackrel{\·}{m}=-\frac{T_{max}u}{g_0{I}_{sp}}\end{array}---\left(1\right)$

其中r＝[x,y,z]为航天器的位置矢量，v＝[v_x,v_y,v_z]为航天器的速度矢量，m为航天器质量，α＝[α_x,α_y,α_z]为发动机推力方向单位矢量，u∈[0,1]为发动机的开关控制量，T_max为发动机最大推力，μ为地球引力常数，g₀为地球水平面重力加速度，I_sp为发动机比冲。

考虑非球形J2摄动对近地卫星运动的影响，由于引力项在三个方向不具有统一性，此时动力学方程表示为公式(2)：

$\begin{array}{c}\stackrel{\·}{x}=v_x\\ \stackrel{\·}{y}=v_y\\ \stackrel{\·}{z}=v_z\\ {\stackrel{\·}{v}}_x=-\frac{\mathrm{\μ}}{r^3}x+\frac{T_{\max }u}{m}{\mathrm{\α}}_x+{\mathrm{\μxJ}}_2\frac{R_E^2}{r^5}(-\frac{3}{2}+\frac{15}{2}\frac{z^2}{r^2})\\ {\stackrel{\·}{v}}_y=-\frac{\mathrm{\μ}}{r^3}y+\frac{T_{\max }u}{m}{\mathrm{\α}}_y+{\mathrm{\μyJ}}_2\frac{R_E^2}{r^5}(-\frac{3}{2}+\frac{15}{2}\frac{z^2}{r^2})\\ {\stackrel{\·}{v}}_z=-\frac{\mathrm{\μ}}{r^3}z+\frac{T_{\max }u}{m}{\mathrm{\α}}_z+{\mathrm{\μzJ}}_2\frac{R_E^2}{r^5}(-\frac{9}{2}+\frac{15}{2}\frac{z^2}{r^2})\\ \stackrel{\·}{m}=-\frac{T_{\max }u}{g_0{I}_{sp}}\end{array}---\left(2\right)$

其中R_E＝6378.137km为地球半径，J₂＝1.082626×10^-3为摄动系数。

步骤二：根据始末航天器状态，求解燃耗最优性能指标所对应的哈密顿函数H。

给定起始时刻t₀航天器的位置矢量[r₀；v₀]，到达时刻t_f的位置矢量为[r_f；v_f]。以燃耗最优作为性能指标J，表达式为公式(3)：

$J=\frac{T_{max}}{g_0{I}_{sp}}{\∫}_{t_0}^{t_f}udt---\left(3\right)$

如果性能指标J中控制变量u为一次的，最优控制率应为Bang-Bang控制。控制变量u应取1或0，控制变量u是不连续的。为使控制变量u连续，通过定义一个过渡系数ε来使控制变量u连续，则构造性能指标J：

$J=\frac{T_{max}}{g_0{I}_{sp}}{\∫}_{t_0}^{t_f}\{u-\mathrm{\ϵ}ln\[u(1-u)\]\}dt---\left(4\right)$

对构造公式(4)的性能指标J乘以一个正的参数λ₀，性能指标J函数形式为公式(5)：

$J={\mathrm{\λ}}_0\frac{T_{\max }}{g_0{I}_{sp}}{\∫}_{t_0}^{t_f}\{u-\mathrm{\ϵ}ln\[u(1-u)\]\}dt---\left(5\right)$

对于构造公式(5)的性能指标J，哈密顿函数H为：

$H={\mathrm{\λ}}_r\·\stackrel{\·}{r}+{\mathrm{\λ}}_v\·\stackrel{\·}{v}-{\mathrm{\λ}}_m\stackrel{\·}{m}+{\mathrm{\λ}}_0\frac{T_{max}}{g_0{I}_{sp}}\{u-\mathrm{\ϵ}ln\[u(1-u)\]\}---\left(6\right)$

步骤三：根据庞特里亚金极大值原理，求解最优控制率以及相应的协态方程。

由庞特里亚金极大值原理可知，为求解最优控制率，需令哈密顿函数最小化，即得最优推力矢量的方向和幅值分别如公式(7)和公式(8)所示：

$\mathrm{\α}=-\frac{{\mathrm{\λ}}_v}{\left|\right|{\mathrm{\λ}}_v\left|\right|}---\left(7\right)$

$u=\frac{2\mathrm{\ϵ}}{\mathrm{\ρ}+2\mathrm{\ϵ}+\sqrt{{\mathrm{\ρ}}^2+4{\mathrm{\ϵ}}^2}}---\left(8\right)$

开关函数ρ为：

$\mathrm{\ρ}=1-\frac{g_0{I}_{sp}\left|\right|{\mathrm{\λ}}_v\left|\right|}{{\mathrm{\λ}}_{0}m}-\frac{{\mathrm{\λ}}_m}{{\mathrm{\λ}}_0}---\left(9\right)$

对于构造公式(5)的性能指标J，虽然最优控制力总是连续并且可导的，不过实际计算中发现过渡系数ε趋于0时，控制力变化仍很剧烈，对协态方程和状态方程积分时精度仍难保证。因此在实际计算过程中，过渡系数ε的取值不宜过小，建议采用过渡系数ε逐渐减小的方式进行迭代计算。

对应的协态方程如公式(10)：

$\left\{\begin{array}{c}{\stackrel{\·}{\mathrm{\λ}}}_x=-\frac{\∂H}{\∂x}={\mathrm{\λ}}_{vx}\frac{\mathrm{\μ}}{r^3}-\frac{3\mathrm{\μ}x}{r^5}r\·{\mathrm{\λ}}_v+{\mathrm{\λ}}_{vx}{\mathrm{\μJ}}_2\frac{R_E^2}{r^5}\[\frac{3}{2}-\frac{15}{2}\frac{x^2}{r^2}-\frac{15}{2}\frac{z^2}{r^2}+\frac{105}{2}\frac{x^2{z}^2}{r^4}\]\\ +{\mathrm{\λ}}_{vy}{\mathrm{\μJ}}_2\frac{R_E^2}{r^5}\[-\frac{15}{2}\frac{xy}{r^2}+\frac{105}{2}\frac{{\mathrm{xyz}}^2}{r^4}\]+{\mathrm{\λ}}_{vz}{\mathrm{\μJ}}_2\frac{R_E^2}{r^5}\[-\frac{45}{2}\frac{xz}{r^2}+\frac{105}{2}\frac{{\mathrm{xz}}^3}{r^4}\]\\ {\stackrel{\·}{\mathrm{\λ}}}_y=-\frac{\∂H}{\∂y}={\mathrm{\λ}}_{vy}\frac{\mathrm{\μ}}{r^3}-\frac{3\mathrm{\μ}y}{r^5}r\·{\mathrm{\λ}}_v+{\mathrm{\λ}}_{vx}{\mathrm{\μJ}}_2\frac{R_E^2}{r^5}\[-\frac{15}{2}\frac{xy}{r^2}+\frac{105}{2}\frac{{\mathrm{xyz}}^2}{r^4}\]\\ +{\mathrm{\λ}}_{vy}{\mathrm{\μJ}}_2\frac{R_E^2}{r^5}\[\frac{3}{2}-\frac{15}{2}\frac{y^2}{r^2}-\frac{15}{2}\frac{z^2}{r^2}+\frac{105}{2}\frac{y^2{z}^2}{r^4}\]+{\mathrm{\λ}}_{vz}{\mathrm{\μJ}}_2\frac{R_E^2}{r^5}\[-\frac{45}{2}\frac{yz}{r^2}+\frac{105}{2}\frac{{\mathrm{yz}}^3}{r^4}\]\\ {\stackrel{\·}{\mathrm{\λ}}}_z=-\frac{\∂H}{\∂z}={\mathrm{\λ}}_{vz}\frac{\mathrm{\μ}}{r^3}-\frac{3\mathrm{\μ}z}{r^5}r\·{\mathrm{\λ}}_v+{\mathrm{\λ}}_{vx}{\mathrm{\μJ}}_2\frac{R_E^2}{r^5}\[-\frac{45}{2}\frac{xz}{r^2}+\frac{105}{2}\frac{{\mathrm{xz}}^3}{r^4}\]\\ +{\mathrm{\λ}}_{vy}{\mathrm{\μJ}}_2\frac{R_E^2}{r^5}\[-\frac{45}{2}\frac{yz}{r^2}+\frac{105}{2}\frac{{\mathrm{yz}}^3}{r^4}\]+{\mathrm{\λ}}_{vz}{\mathrm{\μJ}}_2\frac{R_E^2}{r^5}\[\frac{9}{2}-45\frac{z^2}{r^2}+\frac{105}{2}\frac{z^4}{r^4}\]\\ {\stackrel{\·}{\mathrm{\λ}}}_{vx}=-\frac{\∂H}{\∂v_x}=-{\mathrm{\λ}}_x\\ {\stackrel{\·}{\mathrm{\λ}}}_{vy}=-\frac{\∂H}{\∂v_y}=-{\mathrm{\λ}}_y\\ {\stackrel{\·}{\mathrm{\λ}}}_{vz}=-\frac{\∂H}{\∂v_z}=-{\mathrm{\λ}}_z\\ {\stackrel{\·}{\mathrm{\λ}}}_m=-\frac{\∂H}{\∂m}=-\left|\right|{\mathrm{\λ}}_v\left|\right|\frac{T_{\max }u}{m^2}\end{array}\right.---\left(10\right)$

步骤四：根据变分法最优控制条件，构造两点边值问题方程组。

所述的变分法最优控制条件指

根据Bang-Bang控制的横截条件，当状态变量[r；v；m]边值是固定的，则相应的协态变量[λ_r；λ_v；λ_m]是自由的，当状态变量[r；v；m]为自由的，则协态变量[λ_r；λ_v；λ_m]为零。因此终端时刻的自由状态变量质量对应的协态变量λ_m(t_f)＝0。

因此，燃耗最优控制问题求解变为一个包含多个方程的两点边值问题的求解。

Φ(z)＝[r(t_f)-r_f,v(t_f)-v_f,λ_m(t_f)]^T＝0 (11)

公式(11)称为打靶方程组。把协态变量[λ_r；λ_v；λ_m]与数乘因子λ₀共同称为拉格朗日乘子。则此时最优控制问题等价于对拉格朗日乘子的优化问题。对拉格朗日乘子乘以一个正数不影响燃耗最优控制问题的求解，因此为归一化参数，把拉格朗日乘子向量除以其初值的模的大小。定义新的协态变量为：

$\mathrm{\λ}=\frac{\mathrm{\λ}}{\left|\right|\mathrm{\λ}\left(t_0\right)\left|\right|}$

重新定义后的拉格朗日乘子具有归一化的特点，即||λ(t₀)||＝1。。

因此，把归一化条件||λ(t₀)||＝1加进公式(11)的打靶方程组，新的打靶方程组包含8个方程，如公式(12)：

Φ(z)＝[r(t_f)-r_f,v(t_f)-v_f,λ_m(t_f),||λ(t₀)||-1]^T＝0 (12)

在优化时，给定状态变量[r；v；m]与协态变量[λ_r；λ_v；λ_m]的初始猜测值，对公式(2)的动力学方程和公式(10)的协态方程进行积分，得到满足最优控制条件的终端状态参数[r(t_f),v(t_f),λ_m(t_f)]，代入公式(12)即完成两点边值问题方程组构造。

步骤五：根据航天器始末状态[r₀,v₀,r_f,v_f,m]和飞行时间t，求解两点边值问题打靶方程组。

在求解公式(12)的两点边值问题打靶方程组时，初值采用航天器初始状态和随机生成的协态变量初值，通过非线性规划求解公式(12)的两点边值问题方程组。同时采用不同初值打靶计算，在样本充足的条件下统计可行解中燃耗最少的结果。对燃耗最少的结果采用逐渐减小过渡系数ε的方式迭代计算，直到求解公式(12)的两点边值问题打靶方程组满足预设的精度要求。

为获得较多公式(12)的两点边值问题打靶方程组的可行解，在初始打靶计算时，过渡系数ε取0.01～0.1，而在对最优结果迭代计算时，过渡系数ε取到0.001即能够保证最终结果能够较好的符合Bang-Bang控制。

步骤六：通过步骤五中的求解方法对步骤四中构造的两点边值问题打靶方程组进行求解，得到步骤三中公式(7)所示的最优推力矢量方向α和公式(9)所示的发动机开关函数ρ，以所述的推力矢量方向α和发动机开关函数ρ构成的控制率对卫星变轨进行控制，实现提高近地卫星转移过程的实时性和变轨精度，降低变轨转移所需的燃料消耗。

有益效果：

1、本发明公开的一种考虑J2摄动的近地卫星有限推力最优变轨方法，变轨策略中求解得到的开关函数ρ为与时间有关的函数，是一种实时变轨策略，相比未考虑变轨时长的脉冲控制策略能更为合理的反映变轨控制过程。

2、本发明公开的一种考虑J2摄动的近地卫星有限推力最优变轨方法，该方法基于变分法和庞特里亚金极大值原理，结果满足最优控制条件能够极大降低变轨转移所需的燃料消耗。

3、本发明公开的一种考虑J2摄动的近地卫星有限推力最优变轨方法，该方法由于考虑了非球形J2摄动对转移过程的影响，相比不考虑摄动项的二体动力学方程其变轨精度也得到大大提高。

附图说明：

图1本发明方案流程示意图。

图2地心赤道惯性坐标系的示意图。

图3本发明实例中的转移轨道模型示意图。

图4本发明实例中的推力大小变化曲线图。

图5本发明实例中的俯仰角与偏航角变化曲线图。

图6本发明实例中的推力方向单位矢量变化曲线图。

具体实施方式

为了更好地说明本发明的目的和优点，下面通过对一个近地卫星有限推力变轨实例分析，对本发明做出详细解释。

实施例1：

本发明公开的一种考虑J2摄动的近地卫星有限推力最优变轨方法，包括如下步骤：

步骤一：在地心赤道惯性坐标系中建立考虑J2摄动的航天器动力学方程。

首先考虑引入发动机的开关控制量u，在无摄动条件下笛卡尔坐标系中的有限推力航天器动力学方程公式(1)所示，此时笛卡尔坐标系选为地心赤道惯性坐标系，x轴指向春分点方向，z轴与地球自转轴一致，且指向北，y轴与x轴、z轴满足右手定则。

$\begin{array}{c}\stackrel{\·}{r}=v\\ \stackrel{\·}{v}=-\frac{\mathrm{\μ}}{r^3}r+\frac{T_{max}u}{m}\mathrm{\α}\\ \stackrel{\·}{m}=-\frac{T_{max}u}{g_0{I}_{sp}}\end{array}---\left(1\right)$

其中r＝[x,y,z]为航天器的位置矢量，v＝[v_x,v_y,v_z]为航天器的速度矢量，m为航天器质量，α＝[α_x,α_y,α_z]为发动机推力方向单位矢量，u∈[0,1]为发动机的开关控制量，T_max为发动机最大推力，μ为地球引力常数，g₀为地球水平面重力加速度，I_sp为发动机比冲。这里选取发动机最大推力T_max＝50N；发动机比冲I_sp＝311s。

此时考虑引入非球形J2摄动，由于其引力项在三个方向不具有统一性，故动力学方程表示成以下形式：

$\begin{array}{c}\stackrel{\·}{x}=v_x\\ \stackrel{\·}{y}=v_y\\ \stackrel{\·}{z}=v_z\\ {\stackrel{\·}{v}}_x=-\frac{\mathrm{\μ}}{r^3}x+\frac{T_{\max }u}{m}{\mathrm{\α}}_x+{\mathrm{\μxJ}}_2\frac{R_E^2}{r^5}(-\frac{3}{2}+\frac{15}{2}\frac{z^2}{r^2})\\ {\stackrel{\·}{v}}_y=-\frac{\mathrm{\μ}}{r^3}y+\frac{T_{\max }u}{m}{\mathrm{\α}}_y+{\mathrm{\μyJ}}_2\frac{R_E^2}{r^5}(-\frac{3}{2}+\frac{15}{2}\frac{z^2}{r^2})\\ {\stackrel{\·}{v}}_z=-\frac{\mathrm{\μ}}{r^3}z+\frac{T_{\max }u}{m}{\mathrm{\α}}_z+{\mathrm{\μzJ}}_2\frac{R_E^2}{r^5}(-\frac{9}{2}+\frac{15}{2}\frac{z^2}{r^2})\\ \stackrel{\·}{m}=-\frac{T_{\max }u}{g_0{I}_{sp}}\end{array}---\left(2\right)$

其中R_E＝6378.137km为地球半径，J₂＝1.082626×10^-3为摄动系数。

步骤二：根据始末航天器状态[r₀,v₀,r_f,v_f,m₀]，求解燃耗最优性能指标所对应的哈密顿函数。

这里选取航天器初始质量m₀＝200kg；变轨转移时间为6555.4s；起始时刻t₀航天器状态[r₀,v₀]和到达时刻t_f航天器状态[r_f,v_f]参数如表1所示：

表1航天器始末状态参数

初始航天器状态对应的卫星轨道高度600km，轨道倾角60°；末端航天器状态对应的卫星轨道高度800km，轨道倾角65°。

以燃耗最优为性能指标，其表达式为：

$J=\frac{T_{max}}{g_0{I}_{sp}}{\∫}_{t_0}^{t_f}udt---\left(3\right)$

如果性能指标J中控制变量u为一次的，最优控制率应为Bang-Bang控制。控制变量u应取1或0，控制变量u是不连续的。为使控制变量u连续，通过定义一个过渡系数ε来使控制变量u连续，则构造性能指标J：

$J=\frac{T_{max}}{g_0{I}_{sp}}{\∫}_{t_0}^{t_f}\{u-\mathrm{\ϵ}ln\[u(1-u)\]\}dt---\left(4\right)$

对构造公式(4)的性能指标J乘以一个正的参数λ₀，性能指标J函数形式为公式(5)：

$J={\mathrm{\λ}}_0\frac{T_{\max }}{g_0{I}_{sp}}{\∫}_{t_0}^{t_f}\{u-\mathrm{\ϵ}ln\[u(1-u)\]\}dt---\left(5\right)$

对于构造公式(5)的性能指标J，哈密顿函数H为：

$H={\mathrm{\λ}}_r\·\stackrel{\·}{r}+{\mathrm{\λ}}_v\·\stackrel{\·}{v}-{\mathrm{\λ}}_m\stackrel{\·}{m}+{\mathrm{\λ}}_0\frac{T_{max}}{g_0{I}_{sp}}\{u-\mathrm{\ϵ}ln\[u(1-u)\]\}---\left(6\right)$

步骤三：根据庞特里亚金极大值原理，求解最优控制率以及相应的协态方程。

由庞特里亚金极大值原理可知，为求解最优控制率，需令哈密顿函数最小化，即得最优推力矢量的方向和幅值分别如公式(7)和公式(8)所示：

$\mathrm{\α}=-\frac{{\mathrm{\λ}}_v}{\left|\right|{\mathrm{\λ}}_v\left|\right|}---\left(7\right)$

$u=\frac{2\mathrm{\ϵ}}{\mathrm{\ρ}+2\mathrm{\ϵ}+\sqrt{{\mathrm{\ρ}}^2+4{\mathrm{\ϵ}}^2}}---\left(8\right)$

开关函数ρ为：

$\mathrm{\ρ}=1-\frac{g_0{I}_{sp}\left|\right|{\mathrm{\λ}}_v\left|\right|}{{\mathrm{\λ}}_{0}m}-\frac{{\mathrm{\λ}}_m}{{\mathrm{\λ}}_0}---\left(9\right)$

对于构造公式(5)的性能指标J，虽然最优控制力总是连续并且可导的，不过实际计算中发现过渡系数ε趋于0时，控制力变化仍很剧烈，对协态方程和状态方程积分时精度仍难保证。因此在实际计算过程中，过渡系数ε的取值不宜过小，建议采用过渡系数ε逐渐减小的方式进行迭代计算。

对应的协态方程如公式(10)：

$\left\{\begin{array}{c}{\stackrel{\·}{\mathrm{\λ}}}_x=-\frac{\∂H}{\∂x}={\mathrm{\λ}}_{vx}\frac{\mathrm{\μ}}{r^3}-\frac{3\mathrm{\μ}x}{r^5}r\·{\mathrm{\λ}}_v+{\mathrm{\λ}}_{vx}{\mathrm{\μJ}}_2\frac{R_E^2}{r^5}\[\frac{3}{2}-\frac{15}{2}\frac{x^2}{r^2}-\frac{15}{2}\frac{z^2}{r^2}+\frac{105}{2}\frac{x^2{z}^2}{r^4}\]\\ +{\mathrm{\λ}}_{vy}{\mathrm{\μJ}}_2\frac{R_E^2}{r^5}\[-\frac{15}{2}\frac{xy}{r^2}+\frac{105}{2}\frac{{\mathrm{xyz}}^2}{r^4}\]+{\mathrm{\λ}}_{vz}{\mathrm{\μJ}}_2\frac{R_E^2}{r^5}\[-\frac{45}{2}\frac{xz}{r^2}+\frac{105}{2}\frac{{\mathrm{xz}}^3}{r^4}\]\\ {\stackrel{\·}{\mathrm{\λ}}}_y=-\frac{\∂H}{\∂y}={\mathrm{\λ}}_{vy}\frac{\mathrm{\μ}}{r^3}-\frac{3\mathrm{\μ}y}{r^5}r\·{\mathrm{\λ}}_v+{\mathrm{\λ}}_{vx}{\mathrm{\μJ}}_2\frac{R_E^2}{r^5}\[-\frac{15}{2}\frac{xy}{r^2}+\frac{105}{2}\frac{{\mathrm{xyz}}^2}{r^4}\]\\ +{\mathrm{\λ}}_{vy}{\mathrm{\μJ}}_2\frac{R_E^2}{r^5}\[\frac{3}{2}-\frac{15}{2}\frac{y^2}{r^2}-\frac{15}{2}\frac{z^2}{r^2}+\frac{105}{2}\frac{y^2{z}^2}{r^4}\]+{\mathrm{\λ}}_{vz}{\mathrm{\μJ}}_2\frac{R_E^2}{r^5}\[-\frac{45}{2}\frac{yz}{r^2}+\frac{105}{2}\frac{{\mathrm{yz}}^3}{r^4}\]\\ {\stackrel{\·}{\mathrm{\λ}}}_z=-\frac{\∂H}{\∂z}={\mathrm{\λ}}_{vz}\frac{\mathrm{\μ}}{r^3}-\frac{3\mathrm{\μ}z}{r^5}r\·{\mathrm{\λ}}_v+{\mathrm{\λ}}_{vx}{\mathrm{\μJ}}_2\frac{R_E^2}{r^5}\[-\frac{45}{2}\frac{xz}{r^2}+\frac{105}{2}\frac{{\mathrm{xz}}^3}{r^4}\]\\ +{\mathrm{\λ}}_{vy}{\mathrm{\μJ}}_2\frac{R_E^2}{r^5}\[-\frac{45}{2}\frac{yz}{r^2}+\frac{105}{2}\frac{{\mathrm{yz}}^3}{r^4}\]+{\mathrm{\λ}}_{vz}{\mathrm{\μJ}}_2\frac{R_E^2}{r^5}\[\frac{9}{2}-45\frac{z^2}{r^2}+\frac{105}{2}\frac{z^4}{r^4}\]\\ {\stackrel{\·}{\mathrm{\λ}}}_{vx}=-\frac{\∂H}{\∂v_x}=-{\mathrm{\λ}}_x\\ {\stackrel{\·}{\mathrm{\λ}}}_{vy}=-\frac{\∂H}{\∂v_y}=-{\mathrm{\λ}}_y\\ {\stackrel{\·}{\mathrm{\λ}}}_{vz}=-\frac{\∂H}{\∂v_z}=-{\mathrm{\λ}}_z\\ {\stackrel{\·}{\mathrm{\λ}}}_m=-\frac{\∂H}{\∂m}=-\left|\right|{\mathrm{\λ}}_v\left|\right|\frac{T_{\max }u}{m^2}\end{array}\right.---\left(10\right)$

步骤四：根据变分法最优控制条件，构造两点边值问题方程组。

根据Bang-Bang控制的横截条件，当状态变量[r；v；m]边值是固定的，则相应的协态变量[λ_r；λ_v；λ_m]是自由的，当状态变量[r；v；m]为自由的，则协态变量[λ_r；λ_v；λ_m]为零。因此终端时刻的自由状态变量质量对应的协态变量λ_m(t_f)＝0。

因此，燃耗最优控制问题求解变为一个包含多个方程的两点边值问题的求解。

Φ(z)＝[r(t_f)-r_f,v(t_f)-v_f,λ_m(t_f)]^T＝0 (11)

公式(11)称为打靶方程组。把协态变量[λ_r；λ_v；λ_m]与数乘因子λ₀共同称为拉格朗日乘子。则此时最优控制问题等价于对拉格朗日乘子的优化问题。对拉格朗日乘子乘以一个正数不影响燃耗最优控制问题的求解，因此为归一化参数，把拉格朗日乘子向量除以其初值的模的大小。定义新的协态变量为：

$\mathrm{\λ}=\frac{\mathrm{\λ}}{\left|\right|\mathrm{\λ}\left(t_0\right)\left|\right|}$

重新定义后的拉格朗日乘子具有归一化的特点，即||λ(t₀)||＝1。。

因此，把归一化条件||λ(t₀)||＝1加进公式(11)的打靶方程组，新的打靶方程组包含8个方程，如公式(12)：

Φ(z)＝[r(t_f)-r_f,v(t_f)-v_f,λ_m(t_f),||λ(t₀)||-1]^T＝0 (12)

在优化时，给定状态变量[r；v；m]与协态变量[λ_r；λ_v；λ_m]的初始猜测值，对公式(2)的动力学方程和公式(10)的协态方程进行积分，得到满足最优控制条件的终端状态参数[r(t_f),v(t_f),λ_m(t_f)]，代入公式(12)即完成两点边值问题方程组构造。

步骤五：根据航天器始末状态[r₀,v₀,r_f,v_f,m]和飞行时间t，求解两点边值问题打靶方程组。

在求解公式(12)的两点边值问题打靶方程组时，初值采用航天器初始状态和随机生成的协态变量初值，通过非线性规划求解公式(12)的两点边值问题方程组。同时采用不同初值打靶计算，在样本充足的条件下统计可行解中燃耗最少的结果。对燃耗最少的结果采用逐渐减小过渡系数ε的方式迭代计算，直到求解公式(12)的两点边值问题打靶方程组满足预设的精度要求。

为获得较多公式(12)的两点边值问题打靶方程组的可行解，在初始打靶计算时，过渡系数ε取0.01～0.1，而在对最优结果迭代计算时，过渡系数ε取到0.001即能够保证最终结果能够较好的符合Bang-Bang控制。

这里选取迭代计算的位置精度为1m，速度精度为1m/s。通过求解方程(12)可以得到航天器的整个变轨转移控制过程。整个过程消耗燃料43.67kg，实现轨道高度提升了200km，轨道倾角改变了5°。

步骤六：通过步骤五中的求解方法对步骤四中构造的两点边值问题打靶方程组进行求解，得到步骤三中公式(7)所示的最优推力矢量方向α和公式(9)所示的发动机开关函数ρ，以所述的推力矢量方向α和发动机开关函数ρ控制率对卫星变轨进行控制，实现提高近地卫星转移过程的实时性和变轨精度，降低变轨转移所需的燃料消耗。

以上所述的具体描述，对发明的目的、技术方案和有益效果进行了进一步详细说明，所应理解的是，以上所述仅为本发明的具体实施例，用于解释本发明，并不用于限定本发明的保护范围，凡在本发明的精神和原则之内，所做的任何修改、等同替换、改进等，均应包含在本发明的保护范围之内。

Claims (4)

1.一种考虑J2摄动的近地卫星有限推力最优变轨方法，其特征在于：包括如下步骤：

步骤一：在地心赤道惯性坐标系中建立考虑J2摄动的航天器动力学方程；所述的J2摄动为非球形引力摄动；

首先考虑引入发动机的开关控制量u，在无摄动条件下笛卡尔坐标系中的有限推力航天器动力学方程公式(1)所示，此时笛卡尔坐标系选为地心赤道惯性坐标系，x轴指向春分点方向，z轴与地球自转轴一致，且指向北，y轴与x轴、z轴满足右手定则；

$\begin{array}{c}\stackrel{\·}{r}=v\\ \stackrel{\·}{v}=-\frac{\mathrm{\μ}}{r^3}r+\frac{T_{\max }u}{m}\mathrm{\α}\\ \stackrel{\·}{m}=-\frac{T_{\max }u}{g_0{I}_{sp}}\end{array}---\left(1\right)$

其中r＝[x,y,z]为航天器的位置矢量，v＝[v_x,v_y,v_z]为航天器的速度矢量，m为航天器质量，α＝[α_x,α_y,α_z]为发动机推力方向单位矢量，u∈[0,1]为发动机的开关控制量，T_max为发动机最大推力，μ为地球引力常数，g₀为地球水平面重力加速度，I_sp为发动机比冲；

考虑非球形J2摄动对近地卫星运动的影响，由于引力项在三个方向不具有统一性，此时动力学方程表示为公式(2)：

$\begin{array}{c}\stackrel{\·}{x}=v_x\\ \stackrel{\·}{y}=v_y\\ \stackrel{\·}{z}=v_z\\ {\stackrel{\·}{v}}_x=-\frac{\mathrm{\μ}}{r^3}x+\frac{T_{\max }u}{m}{\mathrm{\α}}_x+{\mathrm{\μxJ}}_2\frac{R_E^2}{r^5}\left(-\frac{3}{2}+\frac{15}{2}\frac{z^2}{r^2}\right)\\ {\stackrel{\·}{v}}_y=-\frac{\mathrm{\μ}}{r^3}y+\frac{T_{\max }u}{m}{\mathrm{\α}}_y+{\mathrm{\μyJ}}_2\frac{R_E^2}{r^5}\left(-\frac{3}{2}+\frac{15}{2}\frac{z^2}{r^2}\right)\\ {\stackrel{\·}{v}}_z=-\frac{\mathrm{\μ}}{r^3}z+\frac{T_{\max }u}{m}{\mathrm{\α}}_z+{\mathrm{\μzJ}}_2\frac{R_E^2}{r^5}\left(-\frac{9}{2}+\frac{15}{2}\frac{z^2}{r^2}\right)\\ \stackrel{\·}{m}=-\frac{T_{\max }u}{g_0{I}_{sp}}\end{array}---\left(2\right)$

其中R_E＝6378.137km为地球半径，J₂＝1.082626×10^-3为摄动系数；

步骤二：根据始末航天器状态，求解燃耗最优性能指标所对应的哈密顿函数H；

给定起始时刻t₀航天器的位置矢量[r₀；v₀]，到达时刻t_f的位置矢量为[r_f；v_f]；以燃耗最优作为性能指标J，表达式为公式(3)：

$J=\frac{T_{max}}{g_0{I}_{sp}}{\∫}_{t_0}^{t_f}udt---\left(3\right)$

如果性能指标J中控制变量u为一次的，最优控制率应为Bang-Bang控制；控制变量u应取1或0，控制变量u是不连续的；为使控制变量u连续，通过定义过渡系数ε来使控制变量u连续，则构造性能指标J：

$J=\frac{T_{max}}{g_0{I}_{sp}}{\∫}_{t_0}^{t_f}\{u-\mathrm{\ϵ}ln\[u(1-u)\]\}dt---\left(4\right)$

对构造公式(4)的性能指标J乘以一个正的参数λ₀，性能指标J函数形式为公式(5)：

$J={\mathrm{\λ}}_0\frac{T_{\max }}{g_0{I}_{sp}}{\∫}_{t_0}^{t_f}\{u-\mathrm{\ϵ}ln\[u(1-u)\]\}dt---\left(5\right)$

对于构造公式(5)的性能指标J，哈密顿函数H为：

$H={\mathrm{\λ}}_r\·\stackrel{\·}{r}+{\mathrm{\λ}}_v\·\stackrel{\·}{v}-{\mathrm{\λ}}_m\stackrel{\·}{m}+{\mathrm{\λ}}_0\frac{T_{max}}{g_0{I}_{sp}}\{u-\mathrm{\ϵ}ln\[u(1-u)\]\}---\left(6\right)$

步骤三：根据庞特里亚金极大值原理，求解最优控制率以及相应的协态方程；

由庞特里亚金极大值原理可知，为求解最优控制率，需令哈密顿函数最小化，即得最优推力矢量的方向和幅值分别如公式(7)和公式(8)所示：

$\mathrm{\α}=-\frac{{\mathrm{\λ}}_v}{\left|\right|{\mathrm{\λ}}_v\left|\right|}---\left(7\right)$

$u=\frac{2\mathrm{\ϵ}}{\mathrm{\ρ}+2\mathrm{\ϵ}+\sqrt{{\mathrm{\ρ}}^2+4{\mathrm{\ϵ}}^2}}---\left(8\right)$

开关函数ρ为：

$\mathrm{\ρ}=1-\frac{g_0{I}_{sp}\left|\right|{\mathrm{\λ}}_v\left|\right|}{{\mathrm{\λ}}_{0}m}-\frac{{\mathrm{\λ}}_m}{{\mathrm{\λ}}_0}---\left(9\right)$

对应的协态方程如公式(10)：

$\left\{\begin{array}{c}{\stackrel{\·}{\mathrm{\λ}}}_x=-\frac{\∂H}{\∂x}={\mathrm{\λ}}_{vx}\frac{\mathrm{\μ}}{r^3}-\frac{3\mathrm{\μ}x}{r^5}r\·{\mathrm{\λ}}_v+{\mathrm{\λ}}_{vx}{\mathrm{\μJ}}_2\frac{R_E^2}{r^5}\[\frac{3}{2}-\frac{15}{2}\frac{x^2}{r^2}-\frac{15}{2}\frac{z^2}{r^2}+\frac{105}{2}\frac{x^2{z}^2}{r^4}\]\\ +{\mathrm{\λ}}_{vy}{\mathrm{\μJ}}_2\frac{R_E^2}{r^5}\[-\frac{15}{2}\frac{xy}{r^2}+\frac{105}{2}\frac{{\mathrm{xyz}}^2}{r^4}\]+{\mathrm{\λ}}_{vz}{\mathrm{\μJ}}_2\frac{R_E^2}{r^5}\[-\frac{45}{2}\frac{xz}{r^2}+\frac{105}{2}\frac{{\mathrm{xz}}^3}{r^4}\]\\ {\stackrel{\·}{\mathrm{\λ}}}_y=-\frac{\∂H}{\∂y}={\mathrm{\λ}}_{vy}\frac{\mathrm{\μ}}{r^3}-\frac{3\mathrm{\μ}y}{r^5}r\·{\mathrm{\λ}}_v+{\mathrm{\λ}}_{vx}{\mathrm{\μJ}}_2\frac{R_E^2}{r^5}\[-\frac{15}{2}\frac{xy}{r^2}+\frac{105}{2}\frac{{\mathrm{xyz}}^2}{r^4}\]\\ +{\mathrm{\λ}}_{vy}{\mathrm{\μJ}}_2\frac{R_E^2}{r^5}\[\frac{3}{2}-\frac{15}{2}\frac{y^2}{r^2}-\frac{15}{2}\frac{z^2}{r^2}+\frac{105}{2}\frac{y^2{z}^2}{r^4}\]+{\mathrm{\λ}}_{vz}{\mathrm{\μJ}}_2\frac{R_E^2}{r^5}\[-\frac{45}{2}\frac{yz}{r^2}+\frac{105}{2}\frac{{\mathrm{yz}}^3}{r^4}\]\\ {\stackrel{\·}{\mathrm{\λ}}}_z=-\frac{\∂H}{\∂z}={\mathrm{\λ}}_{vz}\frac{\mathrm{\μ}}{r^3}-\frac{3\mathrm{\μ}z}{r^5}r\·{\mathrm{\λ}}_v+{\mathrm{\λ}}_{vx}{\mathrm{\μJ}}_2\frac{R_E^2}{r^5}\[-\frac{45}{2}\frac{xz}{r^2}+\frac{105}{2}\frac{{\mathrm{xz}}^3}{r^4}\]\\ +{\mathrm{\λ}}_{vy}{\mathrm{\μJ}}_2\frac{R_E^2}{r^5}\[-\frac{45}{2}\frac{yz}{r^2}+\frac{105}{2}\frac{{\mathrm{yz}}^3}{r^4}\]+{\mathrm{\λ}}_{vz}{\mathrm{\μJ}}_2\frac{R_E^2}{r^5}\[\frac{9}{2}-45\frac{z^2}{r^2}+\frac{105}{2}\frac{z^4}{r^4}\]\\ {\stackrel{\·}{\mathrm{\λ}}}_{vx}=-\frac{\∂H}{\∂v_x}=-{\mathrm{\λ}}_x\\ {\stackrel{\·}{\mathrm{\λ}}}_{vy}=-\frac{\∂H}{\∂v_y}=-{\mathrm{\λ}}_y\\ {\stackrel{\·}{\mathrm{\λ}}}_{vz}=-\frac{\∂H}{\∂v_z}=-{\mathrm{\λ}}_z\\ {\stackrel{\·}{\mathrm{\λ}}}_m=-\frac{\∂H}{\∂m}=-\left|\right|{\mathrm{\λ}}_v\left|\right|\frac{T_{\max }u}{m^2}\end{array}\right.---\left(10\right)$

步骤四：根据变分法最优控制条件，构造两点边值问题方程组；

所述的变分法最优控制条件指

根据Bang-Bang控制的横截条件，当状态变量[r；v；m]边值是固定的，则相应的协态变量[λ_r；λ_v；λ_m]是自由的，当状态变量[r；v；m]为自由的，则协态变量[λ_r；λ_v；λ_m]为零；因此终端时刻的自由状态变量质量对应的协态变量λ_m(t_f)＝0；

因此，燃耗最优控制问题求解变为包含多个方程的两点边值问题的求解；

Φ(z)＝[r(t_f)-r_f，v(t_f)-v_f，λ_m(t_f)]^T＝0 (11)

公式(11)称为打靶方程组；把协态变量[λ_r；λ_v；λ_m]与数乘因子λ₀共同称为拉格朗日乘子；则此时最优控制问题等价于对拉格朗日乘子的优化问题；对拉格朗日乘子乘以一个正数不影响燃耗最优控制问题的求解，因此为归一化参数，把拉格朗日乘子向量除以其初值的模的大小；定义新的协态变量为：

$\mathrm{\λ}=\frac{\mathrm{\λ}}{\left|\right|\mathrm{\λ}\left(t_0\right)\left|\right|}$

重新定义后的拉格朗日乘子具有归一化的特点，即||λ(t₀)||＝1；；

因此，把归一化条件||λ(t₀)||＝1加进公式(11)的打靶方程组，新的打靶方程组包含8个方程，如公式(12)：

Φ(z)＝[r(t_f)-r_f,v(t_f)-v_f,λ_m(t_f),||λ(t₀)||-1]^T＝0 (12)

在优化时，给定状态变量[r；v；m]与协态变量[λ_r；λ_v；λ_m]的初始猜测值，对公式(2)的动力学方程和公式(10)的协态方程进行积分，得到满足最优控制条件的终端状态参数[r(t_f),v(t_f),λ_m(t_f)]，代入公式(12)即完成两点边值问题方程组构造；

步骤五：根据航天器始末状态[r₀,v₀,r_f,v_f,m]和飞行时间t，求解两点边值问题打靶方程组；

在求解公式(12)的两点边值问题打靶方程组时，初值采用航天器初始状态和随机生成的协态变量初值，通过非线性规划求解公式(12)的两点边值问题方程组；同时采用不同初值打靶计算，在样本充足的条件下统计可行解中燃耗最少的结果；对燃耗最少的结果采用逐渐减小过渡系数ε的方式迭代计算，直到求解公式(12)的两点边值问题打靶方程组满足预设的精度要求；

步骤六：通过步骤五中的求解方法对步骤四中构造的两点边值问题打靶方程组进行求解，得到步骤三中公式(7)所示的最优推力矢量方向α和公式(9)所示的发动机开关函数ρ，以所述的推力矢量方向α和发动机开关函数ρ构成的控制率对卫星变轨进行控制，实现提高近地卫星转移过程的实时性和变轨精度，降低变轨转移所需的燃料消耗。

2.如权利要求1所述的一种考虑J2摄动的近地卫星有限推力最优变轨方法，其特征在于：步骤三中对于构造公式(5)的性能指标J，虽然最优控制力总是连续并且可导的，不过实际计算中发现过渡系数ε趋于0时，控制力变化仍很剧烈，对协态方程和状态方程积分时精度仍难保证；因此在实际计算过程中，过渡系数ε的取值不宜过小，采用过渡系数ε逐渐减小的方式进行迭代计算。

3.如权利要求2所述的一种考虑J2摄动的近地卫星有限推力最优变轨方法，其特征在于：为获得较多公式(12)的两点边值问题打靶方程组的可行解，在初始打靶计算时，过渡系数ε取0.01～0.1，而在对最优结果迭代计算时，过渡系数ε取到0.001即能够保证最终结果能够较好的符合Bang-Bang控制。

4.一种考虑J2摄动的近地卫星有限推力最优变轨方法，其特征在于：在地心赤道惯性坐标系中建立考虑J2摄动的航天器动力学方程；根据始末航天器状态，求解燃耗最优性能指标所对应的哈密顿函数H；根据庞特里亚金极大值原理求解最优控制率以及相应的协态方程；根据变分法最优控制条件构造满足燃耗最优的两点边值问题打靶方程组；根据航天器始末状态和飞行时间，对两点边值问题打靶方程组进行求解，得到由最优推力控制方向和发动机开关函数构成的最优控制率；利用所述的控制率对卫星变轨进行控制，实现提高近地卫星转移过程的实时性和变轨精度，降低变轨转移所需的燃料消耗。