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

摘要

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

主权项

一种考虑J2摄动的近地卫星有限推力最优变轨方法，其特征在于：包括如下步骤：步骤一：在地心赤道惯性坐标系中建立考虑J2摄动的航天器动力学方程；所述的J2摄动为非球形引力摄动；首先考虑引入发动机的开关控制量u，在无摄动条件下笛卡尔坐标系中的有限推力航天器动力学方程公式(1)所示，此时笛卡尔坐标系选为地心赤道惯性坐标系，x轴指向春分点方向，z轴与地球自转轴一致，且指向北，y轴与x轴、z轴满足右手定则；$\begin{array}{c}\stackrel{·}{r}=v\\ \stackrel{·}{v}=-\frac{\mu }{r^3}r+\frac{T_{\max }u}{m}\alpha \\ \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{\mu }{r^3}x+\frac{T_{\max }u}{m}{\alpha }_x+{\mathrm{\mu 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{\mu }{r^3}y+\frac{T_{\max }u}{m}{\alpha }_y+{\mathrm{\mu 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{\mu }{r^3}z+\frac{T_{\max }u}{m}{\alpha }_z+{\mathrm{\mu 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}}{\int }_{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}}{\int }_{t_0}^{t_f}\{u-ϵln[u(1-u)\left]\right\}dt---\left(4\right)$对构造公式(4)的性能指标J乘以一个正的参数λ₀，性能指标J函数形式为公式(5)：$J={\lambda }_0\frac{T_{\max }}{g_0{I}_{sp}}{\int }_{t_0}^{t_f}\{u-ϵln[u(1-u)\left]\right\}dt---\left(5\right)$对于构造公式(5)的性能指标J，哈密顿函数H为：$H={\lambda }_r·\stackrel{·}{r}+{\lambda }_v·\stackrel{·}{v}-{\lambda }_m\stackrel{·}{m}+{\lambda }_0\frac{T_{max}}{g_0{I}_{sp}}\{u-ϵln[u(1-u)\left]\right\}---\left(6\right)$步骤三：根据庞特里亚金极大值原理，求解最优控制率以及相应的协态方程；由庞特里亚金极大值原理可知，为求解最优控制率，需令哈密顿函数最小化，即得最优推力矢量的方向和幅值分别如公式(7)和公式(8)所示：$\alpha =-\frac{{\lambda }_v}{\left|\right|{\lambda }_v\left|\right|}---\left(7\right)$$u=\frac{2ϵ}{\rho +2ϵ+\sqrt{{\rho }^2+4{ϵ}^2}}---\left(8\right)$开关函数ρ为：$\rho =1-\frac{g_0{I}_{sp}\left|\right|{\lambda }_v\left|\right|}{{\lambda }_{0}m}-\frac{{\lambda }_m}{{\lambda }_0}---\left(9\right)$对应的协态方程如公式(10)：$\left\{\begin{array}{c}{\stackrel{·}{\lambda }}_x=-\frac{\partial H}{\partial x}={\lambda }_{vx}\frac{\mu }{r^3}-\frac{3\mu x}{r^5}r·{\lambda }_v+{\lambda }_{vx}{\mathrm{\mu J}}_2\frac{R_E^2}{r^5}\left[\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}\right]\\ +{\lambda }_{vy}{\mathrm{\mu J}}_2\frac{R_E^2}{r^5}\left[-\frac{15}{2}\frac{xy}{r^2}+\frac{105}{2}\frac{{\mathrm{xyz}}^2}{r^4}\right]+{\lambda }_{vz}{\mathrm{\mu J}}_2\frac{R_E^2}{r^5}\left[-\frac{45}{2}\frac{xz}{r^2}+\frac{105}{2}\frac{{\mathrm{xz}}^3}{r^4}\right]\\ {\stackrel{·}{\lambda }}_y=-\frac{\partial H}{\partial y}={\lambda }_{vy}\frac{\mu }{r^3}-\frac{3\mu y}{r^5}r·{\lambda }_v+{\lambda }_{vx}{\mathrm{\mu J}}_2\frac{R_E^2}{r^5}\left[-\frac{15}{2}\frac{xy}{r^2}+\frac{105}{2}\frac{{\mathrm{xyz}}^2}{r^4}\right]\\ +{\lambda }_{vy}{\mathrm{\mu J}}_2\frac{R_E^2}{r^5}\left[\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}\right]+{\lambda }_{vz}{\mathrm{\mu J}}_2\frac{R_E^2}{r^5}\left[-\frac{45}{2}\frac{yz}{r^2}+\frac{105}{2}\frac{{\mathrm{yz}}^3}{r^4}\right]\\ {\stackrel{·}{\lambda }}_z=-\frac{\partial H}{\partial z}={\lambda }_{vz}\frac{\mu }{r^3}-\frac{3\mu z}{r^5}r·{\lambda }_v+{\lambda }_{vx}{\mathrm{\mu J}}_2\frac{R_E^2}{r^5}\left[-\frac{45}{2}\frac{xz}{r^2}+\frac{105}{2}\frac{{\mathrm{xz}}^3}{r^4}\right]\\ +{\lambda }_{vy}{\mathrm{\mu J}}_2\frac{R_E^2}{r^5}\left[-\frac{45}{2}\frac{yz}{r^2}+\frac{105}{2}\frac{{\mathrm{yz}}^3}{r^4}\right]+{\lambda }_{vz}{\mathrm{\mu J}}_2\frac{R_E^2}{r^5}\left[\frac{9}{2}-45\frac{z^2}{r^2}+\frac{105}{2}\frac{z^4}{r^4}\right]\\ {\stackrel{·}{\lambda }}_{vx}=-\frac{\partial H}{\partial v_x}=-{\lambda }_x\\ {\stackrel{·}{\lambda }}_{vy}=-\frac{\partial H}{\partial v_y}=-{\lambda }_y\\ {\stackrel{·}{\lambda }}_{vz}=-\frac{\partial H}{\partial v_z}=-{\lambda }_z\\ {\stackrel{·}{\lambda }}_m=-\frac{\partial H}{\partial m}=-\left|\right|{\lambda }_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]与数乘因子λ₀共同称为拉格朗日乘子；则此时最优控制问题等价于对拉格朗日乘子的优化问题；对拉格朗日乘子乘以一个正数不影响燃耗最优控制问题的求解，因此为归一化参数，把拉格朗日乘子向量除以其初值的模的大小；定义新的协态变量为：$\lambda =\frac{\lambda }{\left|\right|\lambda \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)所示的发动机开关函数ρ，以所述的推力矢量方向α和发动机开关函数ρ构成的控制率对卫星变轨进行控制，实现提高近地卫星转移过程的实时性和变轨精度，降低变轨转移所需的燃料消耗。