CN110466804A - 火箭动力下降着陆过程快速轨迹优化方法 - Google Patents

Abstract

火箭动力下降着陆过程快速轨迹优化方法。本发明公开的燃料最优动力下降着陆快速轨迹规划方法，属于火箭制导领域。本发明实现方法为：对动力下降飞行进行动力学建模并量纲归一化，建立三维无量纲动力学方程；把动力学方程的自变量从时间转化为高度；引入动力下降飞行所需的约束，建立燃料最优的动力下降着陆的最优控制问题；将原始最优控制问题中的非线性动力学方程处理成线性的动力学方程；将一部分的非线性性保留，转化为约束。将非凸约束进行凸化，建立一个凸的最优控制问题；将其在非均匀离散点用四阶Runge‑Kutta方法进行离散，建立二阶锥规划问题；通过迭代求解二阶锥规划问题，直到收敛，即实现动力下降着陆飞行的轨迹规划。本发明有高效的优点，能够提高火箭动力下降段着陆安全性和可靠性。

Description

火箭动力下降着陆过程快速轨迹优化方法

技术领域

本发明属于火箭制导领域，涉及火箭垂直着陆的动力下降段的轨迹规划方法，尤其涉及一种基于凸优化的燃料最优的快速的轨迹规划方法。

背景技术

在最近几年中，可重复使用火箭在世界范围内引起了极大地关注。这可以极大地降低火箭发射成本，和更加容易和迅速地执行一次发射任务。在实现火箭的精确着陆过程中,动力下降段制导扮演着极其重要的角色。

燃料最优的动力下降着陆轨迹规划问题是一个典型的最优控制问题。为了提高火箭在动力下降段的抗干扰能力和机动性，需要实时地去求解这个非凸的最优控制问题，以实现火箭的精确着陆。直接通过现有方法(非线性规划或者基于最优控制理论的打靶法)求解这个非凸问题往往是不能应用于实际工程中，因为此类方法无法保证其可靠性，并且求解效率不高。目前对于动力下降着陆的轨迹规划问题的求解，大多数是依赖于直接法求解，例如非线性规划。但这些算法也没法保证求解的可靠性和效率，目前还不能进行工程应用。

发明内容

针对现有的动力下降制导方法的可靠性和效率的不足，本发明公开的一种燃料最优的动力下降段着陆的快速的轨迹规划方法要解决的技术问题是：实现燃料最优的动力下降着陆轨迹的在线规划，具有高效的优点。本发明能够提高火箭动力下降段着陆安全性和可靠性。

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

本发明公开的燃料最优动力下降着陆快速轨迹规划方法，对动力下降飞行进行动力学建模并量纲归一化，建立三维无量纲动力学方程。由于飞行时间是自由的，并把动力学方程的自变量从时间转化为高度。并引入动力下降飞行所需的约束，建立燃料最优的动力下降着陆的最优控制问题。对原始最优控制问题中的非线性动力学方程进行处理，使之成为一个线性的动力学方程。将一部分的非线性性保留，转化为约束。将非凸的约束进行凸化，并建立一个凸的最优控制问题。并将其在非均匀离散点用四阶Runge-Kutta方法进行离散，建立二阶锥规划问题。最后通过迭代求解二阶锥规划问题，直到收敛，即实现动力下降着陆飞行的轨迹规划。

本发明公开的燃料最优快速火箭动力下降着陆轨迹规划方法，包括如下步骤：

步骤一：对动力下降飞行进行动力学建模并量纲归一化，建立三维无量纲动力学方程。由于飞行时间是自由的，并把动力学方程的自变量从时间转化为高度。并引入动力下降飞行所需的约束，建立燃料最优的动力下降着陆的最优控制问题。

步骤一实现方法为：

对火箭动力下降段进行动力学建模，并量纲归一化，火箭动力下降段的无量纲动力学方程表示为

其中，r＝[r₁，r₂，r₃]^T是火箭的空间位置，e₁轴指向高度方向，e₂轴指向正东方向，e₃轴与e₁，e₂构成右手法则；v＝[v₁，v₂，v₃]^T是火箭的速度矢量；m是火箭的质量；g＝[g，0，0]^T是重力加速度矢量，其中g是重力加速度，被考虑为常数；T＝[T₁，T₂，T₃]^T是推力矢量；D表示气动阻力矢量；I_sp是火箭发动机的比冲。在式(1)中，除了角度变量以外，其他变量均进行了量纲归一化，位置变量r用位置因子r_scale归一化，速度v用速度因子v_scale归一化，质量m用质量因子m_scale归一化，时间和比冲I_sp用r_scale/v_scale归一化，重力加速度g用归一化，推力T用归一化。其中无量纲的阻力表示为

其中，ρ是无量纲的空气密度，随高度而变化，S_ref是火箭的无量纲的参考面积，C_D是阻力系数。

火箭在动力下降段的飞行时间是未知的，但是飞行高度是已知的，且在动力下降段，高度是随时间单调下降的，即火箭不会往上飞。因此将式(1)的自变量转化为高度，得

其中k_D＝0.5ρS_refC_D，上标(′)表示对高度r₁求导。用τ代替r₁表示自变量。

接下来，引入动力下降飞行中所需要满足的约束。在动力下降段飞行中，由于火箭的发动机的性能，需要对推力的大小进行如下约束

T_min≤||T||≤T_max (4)

其中，T_min＞0和T_max是可允许的最小最大推力。此外，火箭的倾斜角，即火箭纵轴与e₁轴的夹角，不能大于可允许的最大倾斜角θ_max。此约束等价于对推力方向进行约束，即

||T||cosθ_max≤T₁ (5)

其中0≤θ_max≤90°。另一个重要的约束是飞行斜度约束，即要求火箭飞行在一个锥内。可允许最大的飞行斜度角定义为γ_max，且满足0≤γ_max≤90°。飞行斜度约束被描述为

||r||cosγ_max≤τ (6)

在约束(5)-(6)中，为了使火箭垂直着陆，θ_max和γ_max是随高度变化的。且在靠近着陆点的时候θ_max和γ_max的取值很小，这样可以垂直着陆且避免火箭翻倒。

飞行需要满足的末端条件是

其中τ_f是末端高度，|v_f|是可允许的最大着陆速度。

该问题的优化目标为使燃料消耗最小，所以燃料最优的动力下降着陆的最优控制问题是

P0：min -m(τ_f) (8)

s.t.x′＝f(x，T) (9)

T_min≤||T||≤T_max (10)

||T||cosθ_max≤T₁ (11)

||r||cosγ_max≤τ (12)

Φ(x(τ_f))≤0 (13)

其中，式(9)代表动力学方程(3)，式(13)代表末端约束(7)。在问题P0中，约束(11)-(12)是凸的，但是动力学方程(9)是非线性的，约束(10)是非凸的。因此问题P0是一个非凸问题。

步骤二：对问题P0中的非线性动力学方程进行处理，使之成为一个线性的动力学方程。将动力学方程中一部分非线性性保留，并转化为约束。基于线性的动力学方程，建立新的非凸最优控制问题P1。

步骤二实现方法为：

首先重新定义控制变量如下

通过式(14)-(15)的定义，动力学方程(3)有关推力的项变为线性项。此外，新定义的控制量需满足如下约束

式(16)是因重新定义控制量而额外产生的约束，而且该约束是非凸的。式(16)即是从动力学方程(3)中所保留下来的非线性性。

根据式(15)对u₄的定义，原动力学方程(3)中的第六式改写为

式(17)仍然是非线性的。定义一个新的与质量有关的状态量如下

z：＝ln m (18)

相应地，式(17)变换为线性的形式如下

通过对控制量和状态量的新定义，原非线性的动力学方程(3)变换为如下的控制仿射系统

x′＝f(x)+Bu (20)

其中x＝[r₂，r₃，v₁，v₂，v₃，z]^T为状态变量，u＝[u₁，u₂，u₃，u₄]^T是控制变量，并且

接下来，将动力学方程(20)进行部分线性化，即对非线性部分f(x)进行线性化。线性化后的动力学方程为

x′＝A(x^(k))x+Bu+c(x^(k)) (23)

其中x^(k)表示第k次迭代所求出的状态量，c(x^(k))＝f(x^(k))-A{x^(k))x^(k)，并且

在A(x^(k))中，

根据新的控制量的定义(14)-(15)，推力大小约束(4)转化为

推力方向约束(5)转化为

u₄cosθ_max≤u₁ (27)

此外根据新状态量z的定义(18)，目标函数(8)转化为

J＝-z(τ_f) (28)

为了确保对于f(x)线性化的合理性，需施加如下信赖域约束

|x(τ)-x^(k)(τ)|≤δ (29)

其中δ是一个常矢量，上式是一个分量不等式，即

根据以上讨论，建立一个新的燃料最优的动力下降飞行的最优控制问题：

P1：min -z(τ_f) (30)

s.t.x′＝A(x^(k))x+Bu+c(x^(k)) (31)

|x(τ)-x^(k)(τ)|≤δ (32)

-(T_min/v₁)e^-z≤u₄≤-(T_max/v₁)e^-z (34)

u₄cosθ_max≤u₁ (35)

||r||cosγ_max≤τ (36)

Φ(x(τ_f))≤0 (37)

问题P1是一个非凸的问题，其中约束(33)-(34)是非凸的。约束(33)保留了原动力学方程(3)中的非线性成分。

步骤三：将非凸的约束进行凸化，并建立一个凸的燃料最优动力下降飞行的最优控制问题P2。

步骤三的实现方法为：

首先，将约束(34)进行凸化。将约束(34)两边对于进行线性化，即对1/v₁进行线性化得到

那么约束(34)被转换为

其中，

接下来，将非凸约束(33)凸化为如下形式

约束(40)是一个二阶锥约束，称为约束(33)的松弛约束。相比于约束，约束(40)将其可行域扩大。当且仅当约束(40)活跃，约束(33)和约束(40)才是等价的。

将非凸约束凸化后，建立一个凸的燃料最优的动力下降着陆的最优控制问题：

P2：min -z(τ_f) (41)

s.t.x′＝A(x^(k))x+Bu+c(x^(k)) (42)

|x(τ)-x^(k)(τ)|≤δ (43)

u₄cosθ_max≤u₁ (46)

||r||cosγ_max≤τ (47)

Φ(x(τ_f))≤0 (48)在问题P2中，约束(44)能够保证始终是活跃的，即问题P2的最优解{x^*；u^*}，始终满足

步骤四：计算非均匀离散点，将步骤三得到的凸的最优控制问题P2，在该非均匀离散点处离散。其中动力学方程用四阶Runge-Kutta方法进行离散。离散后，建立一个二阶锥规划问题P3。

步骤四的实现方法为：

在飞行末端，火箭的速度非常小，所以如果对自变量高度进行均匀离散，在末段会产生对时间很大的离散间隔，会产生较大的误差。此外，所述情况还不利于步骤五中的序列求解算法的收敛。因此对自变量高度进行非均匀离散，使其在垂直速度较大的地方稀疏而在垂直速度小的地方稠密。高度被离散为{τ₀，...，τ_N}，其中第i个离散点间隔通过下式进行计算

其中上标(.)表示迭代步，v_1，i是在离散点τ_i的垂直速度。然后求解问题P2，产生新的速度剖面，再将离散点更新。但是持续地改变离散点，不利于序列求解算法的收敛，因此当下列条件满足时，停止更新离散点。

其中，是一个常数。

接下来，为了提高求解精度，选用四阶Runge-Kutta法进行对问题P2中的微分方程(42)进行离散。对于微分方程(42)，用四阶Runge-Kutta法离散后的形式为

其中，

其中并且将式(51)整理后，得到：

H_i-1x_i-1+H_ix_i+G_i-1u_i-1+G_iu_i＝b，i＝1，...，N (56)

其中H_i-1，H_i，G_i-1，G_i和b的详细形式可以从动力学方程(42)获得。定义优化变量为 式(56)可以改写为

My＝d (57)

其中M和d的具体形式从式(56)中直接获得。此外，初始状态约束也包含在式(57)中。

将问题P2在离散点{τ₀，...，τ_N}进行离散，得到如下二阶锥规划问题。

P3：min l^Ty (58)

s.t.My＝d (59)

u_4，icosθ_max≤u_1，i (63)

||r_i||cosγ_max≤τ_i (64)

Φ(x_N)≤0 (65)

其中i＝1，...，N。

步骤五：迭代求解步骤四得到的二阶锥规划问题P3。首先选择初始状态剖面x⁽⁰⁾。在每次迭代中，更新离散点并将上次迭代所得轨迹x^(k)，代入问题P2，并将问题离散求解，得到的新的解x^(k+1)。重复此过程，直到前后两次迭代的解一致，即实现动力下降着陆飞行的轨迹规划。

步骤五的实现方法为：

步骤5.1：设置k＝0，选择初始状态剖面x⁽⁰⁾。

步骤5.2：根据式(49)计算离散点间隔即计算离散点。

步骤5.3：在第k+1次迭代时，将求x(k)代入问题P2并基于离散间隔为 的离散点，离散为问题P3并求解，并得到其解为{x^(k+1)；u^(k+1)}。

步骤5.4：判断如下收敛条件是否满足

这是一个分量形式的不等式，其中∈是一个足够小的常数。如果式(66)满足，执行步骤5.6，否则，执行步骤5.5。

步骤5.5：判断是否满足对于离散点更新的停步准则(50)。如果不满足停步准则(50)，执行步骤5.2，否则，执行步骤5.3。

步骤5.6：原问题P0的解为{x^(k+1)；u^(k+1)}。停止迭代，即实现动力下降着陆飞行的轨迹规划。

有益效果：

1.本发明公开的燃料最优动力下降着陆快速轨迹规划方法，对原动力学方程的自变量从时间转换为高度，能够在优化飞行轨迹的同时求出最优飞行时间，而且能够在优化问题中直接加入与高度相关联的过程约束。

2.本发明公开的燃料最优动力下降着陆快速轨迹规划方法，选取燃料最优为目标函数通过对燃料消耗进行优化，能够得到燃料最优的动力下降着陆轨迹，从而能够降低火箭燃料消耗，利于提高火箭的运载能力。

3.本发明公开的燃料最优动力下降着陆快速轨迹规划方法，对非线性的动力学方程进行处理，保留部分的非线性成分，将此非线性成分转至约束中，并将原始最优控制问题转化为序列二阶锥规划问题，通过序列求解二阶锥规划，能够得到满足过程约束且燃料最优的动力下降着陆轨迹，且具有计算量小，计算快速的特点，能够用于火箭箭载计算机进行实时的轨迹规划。

附图说明

图1是本发明实例中可允许最大倾斜角和飞行斜度角图。

图2是本发明的一种燃料最优动力下降着陆快速轨迹规划算法流程图；

图3是本发明步骤一中坐标系，倾斜角和飞行斜度角的定义图；

图4是本发明步骤三原约束和松弛约束所代表的可行域的示意图；

图5是本发明实例中轨迹和推力方向图；

图6是本发明实例中火箭飞行末端的高度航程图；

图7是本发明实例中倾斜角和飞行斜度角变化图；

图8是本发明实例中比推力和速度图；

图9是本发明实例中气动阻力和推力加速度图。

具体实施方式

为了更好的说明本发明的目的和优点，下面结合附图和实例对发明内容做进一步说明。

为了验证方法的可行性，选择一个火箭动力下降着陆任务进行验证。火箭的初始质量为m₀＝23000kg，干质量为m_dry＝15000kg，发动机比冲为I_sp＝360s，气动阻力系数C_D＝1.3，参考面积为S_ref＝9m²，发动机最大推力为T_max＝300kN。推力可调范围是40％～100％T_max.并且设r_scale＝1000m，v_scale＝100m/s，和m_scale＝m₀。选取离散点数为121。此外，信赖域半径设置为

停步准则参数设置为

离散点更新的停步准则参数为

实例中任务的初始和末端状态见表1。

表1实例中任务的初始状态

此外，对于推力方向约束，可允许的最大倾斜角γ_max设置如图1上曲线所示。当高度大于300m时，可允许的最大倾斜角为30度；当高度小于300m时，可允许的最大倾斜角线性地从30度变至0。对于飞行斜度约束，可允许的最大飞行斜度角设置如图1下曲线所示。当高度大于1250m时，可允许的最大飞行斜度角设置为无穷大，即不施加飞行斜度约束；当高度从1250m变化至1000m时，可允许的最大飞行斜度角从20度线性地变至16度；当高度从1000m变化至100m时，可允许的最大飞行斜度角从16度线性地变至12度；当高度小于100m时，可允许的最大飞行斜度角为常数12度。

如图2所示，本实施例公开的一种燃料最优动力下降着陆快速轨迹规划方法，包括如下步骤：

步骤一：对动力下降飞行进行动力学建模并量纲归一化，建立三维无量纲动力学方程。由于飞行时间是自由的，并把动力学方程的自变量从时间转化为高度。并引入动力下降飞行所需的约束，建立燃料最优的动力下降着陆的最优控制问题。

对火箭动力下降段进行动力学建模，并量纲归一化，火箭动力下降段的无量纲动力学方程表示为

其中，r＝[r₁，r₂，r₃]^T是火箭的空间位置，e₁轴指向高度方向，e₂轴指向正东方向，e₃轴与e₁，e₂构成右手法则，坐标系见图3；v＝[v₁，v₂，v₃]^T是火箭的速度矢量；m是火箭的质量；g＝[g，0，0]^T是重力加速度矢量，其中g＝9.81m/s²是重力加速度；T＝[T₁，T₂，T₃]^T是推力矢量；D表示气动阻力矢量；I_sp＝360s是火箭发动机的比冲。在式(4)中，除了角度变量以外，其他变量均进行了量纲归一化，位置变量r用r_scale＝1000m归一化，速度v用v_scale＝100m/s归一化，质量m用m_scale＝m₀＝23000kg归一化，时间和比冲I_sp用r_scale/v_scale＝10s归一化，重力加速度g用归一化，推力T用 归一化。其中无量纲的阻力表示为

其中，ρ是无量纲的空气密度，随高度而变化，S_ref＝9m²/1000²m²＝9×10^-6是火箭的无量纲的参考面积，C_D＝1.3是阻力系数。

火箭在动力下降段的飞行时间是未知的，但是飞行高度是已知的，且在动力下降段，高度是随时间单调下降的，即火箭不会往上飞。因此将式(4)的自变量转化为高度，得

其中k_D＝0.5ρS_refC_D，上标(′)表示对高度r₁求导。用τ代替r₁表示自变量。

接下来，引入动力下降飞行中所需要满足的约束。在动力下降段飞行中，由于火箭的发动机的性能，需要对推力的大小进行如下约束：

T_min≤||T||≤T_max (7)

其中，T_min＝0.4T_max＝120kN和T_max＝300kN是可允许的最小最大推力。此外，如图3所示，火箭的倾斜角，即火箭纵轴与e₁轴的夹角，不能大于可允许的最大倾斜角θ_max。此约束等价于对推力方向进行约束，即

||T||cosθ_max≤T₁ (8)

其中θ_max如图1上所示。另一个重要的约束是飞行斜度约束，即要求火箭飞行在一个锥内，见图3。可允许最大的飞行斜度角定义为γ_max，其值如图1下所示。飞行斜度约束被描述为

||r||cosγ_max≤τ (9)

在约束(8)-(9)中，为了使火箭垂直着陆，θ_max和γ_max是随高度变化的，如图1所示。且在靠近着陆点的时候θ_max和γ_max的取值很小，这样可以垂直着陆且避免火箭翻倒。

飞行需要满足的末端条件是

其中τ_f＝0是末端高度。

该问题的优化目标为使燃料消耗最小，所以燃料最优的动力下降着陆的最优控制问题是

P0：min -m(τ_f) (11)

s.t.x′＝f(x，T) (12)

T_min≤||T||≤T_max (13)

||T||cosθ_max≤T₁ (14)

||r||cosγ_max≤τ (15)

Φ(x(τ_f))≤0 (16)

其中，式(12)代表动力学方程(6)，式(16)代表末端约束(10)。在问题P0中，约束(14)-(15)是凸的，但是动力学方程(12)是非线性的，约束(13)是非凸的。因此问题P0是一个非凸问题。

步骤二：对问题P0中的非线性动力学方程进行处理，使之成为一个线性的动力学方程。将一部分的非线性性保留，转化为约束。基于线性的动力学方程，建立新的非凸最优控制问题P1。

首先重新定义控制变量如下

通过式(17)-(18)的定义，动力学方程(3)有关推力的项变为线性项。此外，新定义的控制量需满足如下约束

式(19)是因重新定义控制量而额外产生的约束，而且该约束是非凸的。式(16)即是从动力学方程(6)中所保留下来的非线性性。

根据式(18)对u₄的定义，原动力学方程(6)中的第六式改写为

式(20)仍然是非线性的。定义一个新的与质量有关的状态量如下，

z：＝ln m (21)

相应地，式(20)变换为线性的形式如下

通过对控制量和状态量的新定义，原非线性的动力学方程(6)变换为如下的控制仿射系统

x′＝f(x)+Bu (23)

其中x＝[r₂，r₃，v₁，v₂，v₃，z]^T为状态变量，u＝[u₁，u₂，u₃，u₄]^T是控制变量，并且

接下来，将动力学方程(23)进行部分线性化，即对非线性部分f(x)进行线性化。线性化后的动力学方程为

x′＝A(x^(k))x+Bu+c(x^(k)) (26)

其中x^(k)表示第k次迭代所求出的状态量，c(x^(k))＝f(x^(k))-A{x^(k))x^(k)，并且

在A(x^(k))中，

根据新的控制量的定义(17)-(18)，推力大小约束(7)转化为

推力方向约束(8)转化为

u₄cosθ_max≤u₁ (30)

此外根据新状态量z的定义(21)，目标函数(11)转化为

J＝-z(τ_f) (31)

为了确保对于f(x)线性化的合理性，需施加如下信赖域约束

|x(τ)-x^(k)(τ)|≤δ (32)

其中δ是一个常矢量，其值见式(1)，上式是一个分量不等式，即

根据以上讨论，建立一个新的燃料最优的动力下降飞行的最优控制问题

P1：min -z(τ_f) (33)

s.t.x′＝A(x^(k))x+Bu+c(x^(k)) (34)

|x(τ)-x^(k)(τ)|≤δ (35)

-(T_min/v₁)e^-z≤u₄≤-(T_max/v₁)e^-z (37)

u₄cosθ_max≤u₁ (38)

||r||cosγ_max≤τ (39)

Φ(x(τ_f))≤0 (40)

问题P1是一个非凸的问题，其中约束(36)-(37)是非凸的。约束(36)保留了原动力学方程(6)中的非线性成分。

步骤三：将非凸的约束进行凸化，并建立一个凸的燃料最优动力下降飞行的最优控制问题P2。

首先，将约束(37)进行凸化。将约束(37)两边对于进行线性化，即对1/v₁进行线性化得到：

那么约束(37)被转换为

其中，

接下来，将非凸约束(36)凸化为如下形式

约束(43)是一个二阶锥约束，称为约束(36)的松弛约束。相比于约束(36)，约束(43)将其可行域扩大。约束(36)和(43)，连同约束(42)，在三维的情况下的可行域如图4所示。当且仅当约束(43)活跃，即求得问题的解位于可行域曲面上，约束(36)和约束(43)才是等价的。

将非凸约束凸化后，建立一个凸的燃料最优的动力下降着陆的最优控制问题

P2：min -z(τ_f) (44)

s.t.x′＝A(x^(k))x+Bu+c(x^(k)) (45)

|x(τ)-x^(k)(τ)|≤δ (46)

u₄cosθ_max≤u₁ (49)

||r||cosγ_max≤τ (50)

Φ(x(τ_f))≤0 (51)

在问题P2中，约束(47)能够保证始终是活跃的，即问题P2的最优解{x^*；u^*}，始终满足

步骤四：计算非均匀离散点，将步骤三得到的凸的最优控制问题P2，在该非均匀离散点处离散。其中动力学方程用四阶Runge-Kutta方法进行离散。离散后，建立一个二阶锥规划问题P3。

在飞行末端，火箭的速度非常小，所以如果对自变量高度进行均匀离散，在末段会产生对时间很大的离散间隔，会产生较大的误差。此外，所述情况还不利于步骤五中的序列求解算法的收敛。因此对自变量高度进行非均匀离散，使其在垂直速度较大的地方稀疏而在垂直速度小的地方稠密。高度被离散为{τ₀，...，τ_N}，其中第i个离散点间隔通过下式进行计算

其中上标(.)表示迭代步，v_1，i是在离散点τ_i的垂直速度。然后求解问题P2，产生新的速度剖面，再将离散点更新。但是持续地改变离散点，不利于序列求解算法的收敛，因此当下列条件满足时，停止更新离散点。

其中，是一个常数，其值见式(3)。

接下来，为了提高求解精度，选用四阶Runge-Kutta法进行对问题P2中的微分方程(45)进行离散。对于微分方程(45)，用四阶Runge-Kutta法离散后的形式为

其中，

其中并且将式(54)整理后，得到：

H_i-1x_i-1+H_ix_i+G_i-1u_i-1+G_iu_i＝b，i＝1，...，N (59)

其中H_i-1，H_i，G_i-1，G_i和b的详细形式可以从动力学方程(45)获得。定义优化变量为 式(59)可以改写为

My＝d (60)

其中M和d的具体形式从式(59)中直接获得。此外，初始状态约束也包含在式(60)中。

将问题P2在离散点{τ₀，...，τ_N}进行离散，得到如下二阶锥规划问题。

P3：min l^Ty (61)

s.t.My＝d (62)

u_4，icosθ_max≤u_1，i (66)

||r_i||cosγ_max≤τ_i (67)

Φ(x_N)≤0 (68)

其中i＝1，...，N。

步骤五：迭代求解步骤四得到的二阶锥规划问题P3。首先选择初始状态剖面x⁽⁰⁾。在每次迭代中，更新离散点并将上次迭代所得轨迹x^(k)，代入问题P2，并将问题离散求解，得到的新的解x^(k+1)。重复此过程，直到前后两次迭代的解一致，即实现动力下降着陆飞行的轨迹规划。

步骤5.1：设置k＝0，选择初始状态剖面x⁽⁰⁾为初始值至末端值的线性函数，其中在这里设置质量的末端值为m_dry＝15000kg。

步骤5.2：根据式(52)计算离散点间隔即计算离散点。

步骤5.3：在第k+1次迭代时，将求x(k)代入问题P2并基于离散间隔为 的离散点，离散为问题P3并求解，并得到其解为{x^(k+1)；u^(k+1)}。

步骤5.4：判断如下收敛条件是否满足

这是一个分量形式的不等式，其中∈是一个足够小的常数。如果式(69)满足，执行步骤5.6，否则，执行步骤5.5。

步骤5.5：判断是否满足对于离散点更新的停步准则(53)。如果不满足停步准则(53)，执行步骤5.2，否则，执行步骤5.3。

步骤5.6：原问题P0的解为{x^(k+1)；u^(k+1)}。停止迭代，即实现动力下降着陆飞行的轨迹规划。

基于该任务，分别求出施加和不施加过程约束(49)-(50)的情况。施加过程约束(49)-(50)的情况，求解过程仅用7次迭代，燃料消耗为3264.7kg，飞行时间为49.1s；不施加过程约束(49)-(50)的情况，求解过程用6次迭代，燃料消耗为2840.9kg，飞行时间为42.5s。施加过程约束(49)-(50)会多消耗14.9％的燃料，增加15.5％的飞行时间。图5展示了这两种情况的三维轨迹图和推力方向图，由此可见，当施加过程约束(49)-(50)的时候，火箭几乎是以垂直的姿态着陆的。图6展示了火箭的末端飞行的高度航程的变化，由此也可见，当有过程约束(49)-(50)时，火箭在末端几乎垂直下降。图7展示了火箭在飞行过程中倾斜角和飞行斜度角的变化规律，当施加过程约束(49)-(50)时，倾斜角和飞行斜度角均在约束范围内。图8上展示了比推力大小图，由此可见，不论过程约束(49)-(50)是否存在，推力大小均为bang-bang结构。图8下展示了飞行速度大小的变化。图9展示了气动阻力加速度和推力加速度的变化。

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

Claims (6)

1.燃料最优动力下降着陆快速轨迹规划方法，其特征在于：包括如下步骤，

步骤一：对动力下降飞行进行动力学建模并量纲归一化，建立三维无量纲动力学方程；由于飞行时间是自由的，并把动力学方程的自变量从时间转化为高度；并引入动力下降飞行所需的约束，建立燃料最优的动力下降着陆的最优控制问题；

步骤二：对问题P0中的非线性动力学方程进行处理，使之成为一个线性的动力学方程；将动力学方程中一部分非线性性保留，并转化为约束；基于线性的动力学方程，建立新的非凸最优控制问题P1；

步骤三：将非凸的约束进行凸化，并建立一个凸的燃料最优动力下降飞行的最优控制问题P2；

步骤四：计算非均匀离散点，将步骤三得到的凸的最优控制问题P2，在该非均匀离散点处离散；其中动力学方程用四阶Runge-Kutta方法进行离散；离散后，建立一个二阶锥规划问题P3；

步骤五：迭代求解步骤四得到的二阶锥规划问题P3；首先选择初始状态剖面x⁽⁰⁾；在每次迭代中，更新离散点并将上次迭代所得轨迹x^(k)，代入问题P2，并将问题离散求解，得到的新的解x^(k+1)；重复此过程，直到前后两次迭代的解一致，即实现动力下降着陆飞行的轨迹规划。

2.如权利要求1所述的燃料最优动力下降着陆快速轨迹规划方法，其特征在于：步骤一实现方法为，

对火箭动力下降段进行动力学建模，并量纲归一化，火箭动力下降段的无量纲动力学方程表示为

其中，r＝[r₁，r₂，r₃]^T是火箭的空间位置，e₁轴指向高度方向，e₂轴指向正东方向，e₃轴与e₁，e₂构成右手法则；v＝[v₁，v₂，v₃]^T是火箭的速度矢量；m是火箭的质量；g＝[g，0，0]^T是重力加速度矢量，其中g是重力加速度，被考虑为常数；T＝[T₁，T₂，T₃]^T是推力矢量；D表示气动阻力矢量；I_sp是火箭发动机的比冲；在式(1)中，除了角度变量以外，其他变量均进行了量纲归一化，位置变量r用位置因子r_scale归一化，速度v用速度因子v_scale归一化，质量m用质量因子m_scale归一化，时间和比冲I_sp用r_scale/v_scale归一化，重力加速度g用归一化，推力T用归一化；其中无量纲的阻力表示为

其中，ρ是无量纲的空气密度，随高度而变化，S_ref是火箭的无量纲的参考面积，C_D是阻力系数；

火箭在动力下降段的飞行时间是未知的，但是飞行高度是已知的，且在动力下降段，高度是随时间单调下降的，即火箭不会往上飞；因此将式(1)的自变量转化为高度，得

其中k_D＝0.5ρS_refC_D，上标(′)表示对高度r₁求导；用τ代替r₁表示自变量；

接下来，引入动力下降飞行中所需要满足的约束；在动力下降段飞行中，由于火箭的发动机的性能，需要对推力的大小进行如下约束

T_min≤||T||≤T_max (4)

其中，T_min＞0和T_max是可允许的最小最大推力；此外，火箭的倾斜角，即火箭纵轴与e₁轴的夹角，不能大于可允许的最大倾斜角θ_max；此约束等价于对推力方向进行约束，即

||T||cosθ_max≤T₁ (5)

其中0≤θ_max≤90°；另一个重要的约束是飞行斜度约束，即要求火箭飞行在一个锥内；可允许最大的飞行斜度角定义为γ_max，且满足0≤γ_max≤90°；飞行斜度约束被描述为

||r||cosγ_max≤τ (6)

在约束(5)-(6)中，为了使火箭垂直着陆，θ_max和γ_max是随高度变化的；且在靠近着陆点的时候θ_max和γ_max的取值很小，这样可以垂直着陆且避免火箭翻倒；

飞行需要满足的末端条件是

其中τ_f是末端高度，|v_f|是可允许的最大着陆速度；

该问题的优化目标为使燃料消耗最小，所以燃料最优的动力下降着陆的最优控制问题是

P0：min -m(τ_f) (8)

s.t. x′＝f(x，T) (9)

T_min≤||T||≤T_max (10)

||T||cosθ_max≤T₁ (11)

||r||cosγ_max≤τ (12)

Φ(x(τ_f))≤0 (13)

其中，式(9)代表动力学方程(3)，式(13)代表末端约束(7)；在问题P0中，约束(11)-(12)是凸的，但是动力学方程(9)是非线性的，约束(10)是非凸的；因此问题P0是一个非凸问题。

3.如权利要求2所述的燃料最优动力下降着陆快速轨迹规划方法，其特征在于：步骤二实现方法为，

首先重新定义控制变量如下

通过式(14)-(15)的定义，动力学方程(3)有关推力的项变为线性项；此外，新定义的控制量需满足如下约束

式(16)是因重新定义控制量而额外产生的约束，而且该约束是非凸的；式(16)即是从动力学方程(3)中所保留下来的非线性性；

根据式(15)对u₄的定义，原动力学方程(3)中的第六式改写为

式(17)仍然是非线性的；定义一个新的与质量有关的状态量如下

z：＝ln m (18)

相应地，式(17)变换为线性的形式如下

通过对控制量和状态量的新定义，原非线性的动力学方程(3)变换为如下的控制仿射系统

x′＝f(x)+Bu (20)

其中x＝[r₂，r₃，v₁，v₂，v₃，z]^T为状态变量，u＝[u₁，u₂，u₃，u₄]^T是控制变量，并且

接下来，将动力学方程(20)进行部分线性化，即对非线性部分f(x)进行线性化；线性化后的动力学方程为

x′＝A(x^(k))x+Bu+c(x^(k)) (23)

其中x^(k)表示第k次迭代所求出的状态量，c(x^(k))＝f(x^(k))-A{x^(k))x^(k)，并且

在A(x^(k))中，

根据新的控制量的定义(14)-(15)，推力大小约束(4)转化为

推力方向约束(5)转化为

u₄cosθ_max≤u₁ (27)

此外根据新状态量z的定义(18)，目标函数(8)转化为

J＝-z(τ_f) (28)

为了确保对于f(x)线性化的合理性，需施加如下信赖域约束

|x(τ)-x^(k)(τ)|≤δ (29)

其中δ是一个常矢量，上式是一个分量不等式，即

根据以上讨论，建立一个新的燃料最优的动力下降飞行的最优控制问题：

P1：min -z(τ_f) (30)

s.t. x′＝A(x^(k))x+Bu+c(x^(k)) (31)

|x(τ)-x^(k)(τ)|≤δ (32)

-(T_min/v₁)e^-z≤u₄≤-(T_max/v₁)e^-z (34)

u₄cosθ_max≤u₁ (35)

||r||cosγ_max≤τ (36)

Φ(x(τ_f))≤0 (37)

问题P1是一个非凸的问题，其中约束(33)-(34)是非凸的；约束(33)保留了原动力学方程(3)中的非线性成分。

4.如权利要求3所述的燃料最优动力下降着陆快速轨迹规划方法，其特征在于：步骤三的实现方法为，

首先，将约束(34)进行凸化；将约束(34)两边对于进行线性化，即对1/v₁进行线性化得到

那么约束(34)被转换为

其中，

接下来，将非凸约束(33)凸化为如下形式

约束(40)是一个二阶锥约束，称为约束(33)的松弛约束；相比于约束，约束(40)将其可行域扩大；当且仅当约束(40)活跃，约束(33)和约束(40)才是等价的；

将非凸约束凸化后，建立一个凸的燃料最优的动力下降着陆的最优控制问题

P2：min -z(τ_f) (41)

s.t. x′＝A(x^(k))x+Bu+c(x^(k)) (42)

|x(τ)-x^(k)(τ)|≤δ (43)

u₄cosθ_max≤u₁ (46)

||r||cosγ_max≤τ (47)

Φ(x(τ_f))≤0 (48)

在问题P2中，约束(44)能够保证始终是活跃的，即问题P2的最优解{x^*；u^*}，始终满足

5.如权利要求4所述的燃料最优快速火箭动力下降着陆轨迹规划方法，其特征在于：步骤四的实现方法为，

在飞行末端，火箭的速度非常小，所以如果对自变量高度进行均匀离散，在末段会产生对时间很大的离散间隔，会产生较大的误差；此外，所述情况还不利于步骤五中的序列求解算法的收敛；因此对自变量高度进行非均匀离散，使其在垂直速度较大的地方稀疏而在垂直速度小的地方稠密；高度被离散为{τ₀，...，τ_N}，其中第i个离散点间隔通过下式进行计算

其中上标(.)表示迭代步，v_1，i是在离散点τ_i的垂直速度；然后求解问题P2，产生新的速度剖面，再将离散点更新；但是持续地改变离散点，不利于序列求解算法的收敛，因此当下列条件满足时，停止更新离散点；

其中，是一个常数；

接下来，为了提高求解精度，选用四阶Runge-Kutta法进行对问题P2中的微分方程(42)进行离散；对于微分方程(42)，用四阶Runge-Kutta法离散后的形式为

其中，

其中并且将式(51)整理后，得到

H_i-1x_i-1+H_ix_i+G_i-1u_i-1+G_iu_i＝b，i＝1，...，N (56)

其中H_i-1，H_i，G_i-1，G_i和b的详细形式可以从动力学方程(42)获得；定义优化变量为 式(56)改写为

My＝d (57)

其中M和d的具体形式从式(56)中直接获得；此外，初始状态约束也包含在式(57)中；

将问题P2在离散点{τ₀，...，τ_N}进行离散，得到如下二阶锥规划问题；

P3：min l^Ty (58)

s.t.My＝d (59)

u_4，icosθ_max≤u_1，i (63)

||r_i||cosγ_max≤τ_i (64)

Φ(x_N)≤0 (65)

其中i＝1，...，N。

6.如权利要求5所述的燃料最优快速火箭动力下降着陆轨迹规划方法，其特征在于：步骤五实现方法为，

步骤5.1：设置k＝0，选择初始状态剖面x⁽⁰⁾；

步骤5.2：根据式(49)计算离散点间隔即计算离散点；

步骤5.3：在第k+1次迭代时，将求x^(k)代入问题P2并基于离散间隔为 的离散点，离散为问题P3并求解，并得到其解为{x^(k+1)；u^(k+1)}；

步骤5.4：判断如下收敛条件是否满足

这是一个分量形式的不等式，其中∈是一个足够小的常数；如果式(66)满足，执行步骤5.6，否则，执行步骤5.5；

步骤5.5：判断是否满足对于离散点更新的停步准则(50)；如果不满足停步准则(50)，执行步骤5.2，否则，执行步骤5.3；

步骤5.6：原问题P0的解为{x^(k+1)；u^(k+1)}；停止迭代，即实现动力下降着陆飞行的轨迹规划。