梯度风场中无人机动态滑翔飞行轨迹优化

朱 熠,李继广,郝向宇

(1.陆军工程大学 野战工程学院,南京 210007;2.西安航空学院 飞行器学院,西安 710077)

动态滑翔也称为梯度风场翱翔。一些鸟类,如信天翁,可以长时间利用梯度风场沿着海岸滑翔,它们不会扇动翅膀。自19世纪80年代,学者们开始研究动态飙升的过程,取得了巨大成就,人们关注的焦点是最优轨道的设计和风能的获取。

梯度风场中的动态滑翔对风场的特性有一定的要求。通常,风场的方向要求是水平的,高度方向上是梯度变化。动态滑翔的典型飞行轨迹如图1所示。一个完整的动力翱翔过程一般包括四个阶段,即逆风爬升阶段、高空转弯阶段、顺风潜水阶段和低空转弯阶段。势能在逆风爬升阶段增加,而动能在顺风潜水阶段增加[1]。部分文献表明,风能主要在逆风爬升阶段和高空转弯早期获得,也有一些在顺风潜水阶段获得能量;而部分文献表明,只有逆风爬升和顺风潜水才是风能获取的阶段。

图1 动态滑翔的典型飞行轨迹

在梯度风场中,可以设计无人机的飞行轨迹来优化无人机系统的总能量,考虑到无人机的横向运动,最优飞行轨迹较为复杂[2]。参考文献[3]分析了无人机在无动力飞行中产生最佳轨迹时受不同参数的影响。参考文献[4]对有关动态飙升的文献进行了综述。参考文献[5]研究了高空动态滑翔的技术。

在众多优化算法中,伪谱方法是最近发展起来的一种求解最优控制问题的方法。该方法结合了间接方法和直接方法的优点:它不仅可以获得最优控制问题的数值解,还可以提供协变量的精确信息[6-7]。通过这种方式,连续系统的最优控制问题被转化为受一系列代数约束的非线性规划(Nonlinear programming,NLP)问题。最常用的伪谱方法包括:勒让德(Legendre)伪谱方法、高斯(Gauss)伪谱方法和拉道(Radau)伪谱方法。

本文考虑非常规条件下的动态滑翔轨迹优化过程,采用一种大型无人机进行动态飞行轨迹优化,并用NED坐标系代替常用的ENU坐标系推导了动态腾空的运动方程,用Radau伪谱方法进行轨迹优化。无动力飞行和有动力飞行的最优轨迹都是通过能量优化生成的。

1 动态滑翔的数学描述

1.1 风场模型和惯性坐标系

(1)

因此,为了获得与其他文献相一致的等效表达式,假设风场的方向为y轴方向,即Vwx=Vwz=0.梯度风场通常有三种形式的模型,即对数形式、指数形式和线性形式,Vwy满足以下形式:

(1)对数模型

(2)

其中:Vwy(h)是在高度h时的风速;Vref是参考高度href时的风速;h0是校正系数,代表地面的起伏程度。

(2)指数模型

(3)

(3)线性模型

在公式(3)中, 令p=1可得线性化模型:

(4)

1.2 动态滑翔的运动方程

在本节中考虑了无动力飞行的情况,使用了无人机在飞行轨迹坐标系中的质点模型,纵向和横向运动都被考虑在内。设γ是爬升俯冲角,φ是偏航角,μ是航迹滚转角,它们在飞行轨迹坐标系中。L是无人机的升力,D是无人机的阻力,m是无人机的质量。根据牛顿第二定律和机身坐标系与惯性坐标系之间的转换规则,可得:

(5)

(6)

(7)

因此,无人机在风场中的运动方程可以表示为:

(8)

(9)

(10)

通过将公式(8)至公式(10)分别对时间t求导,可得:

(11)

(12)

(13)

(14)

(15)

(16)

再对公式(14)和公式(15)两边分别乘以cosφ和sinφ, 得到:

(17)

(18)

然后,将公式(17)和公式(18)两边分别相加,可得:

(19)

公式(19)两边同乘以cosγ,可得:

(20)

公式(16)两边同乘以sinγ,得到:

(21)

公式(20)减去公式(21):

(22)

即:

(23)

公式(19)两边同乘以sinγ,公式(16)两边同乘以cosγ, 可得:

(24)

(25)

公式(24)和公式(25)相加:

(26)

即:

(27)

最后,把公式(23)和公式(27)替换到公式(14)和公式(15),可得:

(28)

因为是在NED坐标系中,有h=-zi,因此公式(10)改写为

(29)

公式(8)、(9)、(29)、(23)、(27)、(28)这六个微分方程就构成了状态空间方程,状态变量为[xi,yi,h,Vt,γ′,φ]:

(30)

2 伪谱优化算法的设计

Radau伪谱方法(RPM)求解最优控制问题的基本原理是:将未知状态变量和控制变量离散在一系列Legendre-Gauss-Radau(LGR)点上;构造全局插值多项式来逼近状态变量和控制变量;然后,用状态变量的导数代替动态微分方程[8],具体设计过程如下。

2.1 时域区间转换

伪谱方法最常见的三种点集配置方法是:Legendre-Gauss (LG), Legendre-Gauss-Radau (LGR), Legendre-Gauss-Lobatto (LGL)[9],其中LG、LGR和LGL分别对应于开放区间τ∈(-1,1)、半闭区间τ∈[-1,1)和(-1,1]、闭区间τ∈[-1,1]。无论哪一组点,都需要实现时域区间从t∈[t0,tf]到τ∈[-1,1]。将时域区间t∈[t0,tf]划分成K个网格区间[tk-1,tk],k=1,…,K,其中,t0,…,tK是网格点,并且满足:

t0=t1<…

(31)

在每个网格区间t∈[tk-1,tk]中,令

(32)

由公式(32),区间t∈[tk-1,tk]被转换到τ∈[-1,1],且有:

(33)

2.2 最优控制策略

令x(k)(τ)和u(k)(τ)分别表示状态变量和控制变量的第k个网格点,且τ∈[-1,1]。代价函数的表达式如下:

(34)

具有动态约束、不等式约束和边界条件的方程为:

(35)

C[x(k)(τ),u(k)(τ),τ,tk-1,tk]≤0

(36)

Φ[x(1)(-1),t0,x(K)(1),tf]=0

(37)

2.3 多项式逼近

在处理RPM的连续时间最优控制问题时,有必要根据以下规则使用拉格朗日多项式来近似每个网格区间k∈[1,…,K]上的连续最优问题:

(38)

(39)

(40)

2.4 多项式求导

在用全局插值多项式逼近状态变量之后,状态变量的导数可以用拉格朗日多项式关于参数τ的导数来近似。因此,动态微分方程被转换为代数约束。对公式(38)中的X(k)(τ)关于参数τ求导:

(41)

动态约束方程公式(35)通过公式(41)的Nk个LGR点进行配置:

(42)

(43)

(44)

为了确保网格点的连续性,还需要满足以下等式:

(45)

由代数约束方程公式(42)至公式(45),利用Radau伪谱近似方法,将连续最优控制问题转化为一个最小化代价函数公式(41)的非线性规划问题。

3 仿真分析

3.1 最小风场无动力滑翔

在无动力飞行中,无人机的运动方程如公式(30)。状态变量和控制变量分别为[xi,yi,h,Vt,γ,φ]和[CL,μ]。阻力系数CD没有额外考虑,因为CD和升力系数CL之间的存在如下关系:

(46)

其中,CD0是寄生阻力系数。Rmax是最大升阻比,K是诱导阻力系数,它是CD0和Rmax的函数:

(47)

(48)

表1 样本无人机的参数值

表1中,xi(tf)、yi(tf)和φ(tf)的边界条件并未给出,这是因为通过这三个约束条件的不同设置,可以获得以下三种特殊的飞行状态。

(1)飞行状态1

边界约束:φ(tf)=φ(t0)+π/2,xi(tf)和yi(tf)无约束。飞行轨迹被称为动态翱翔的基本轨迹,航向的变化是90°,φ的范围是[-π/2,π/2]。限制φ范围的目的是设置搜索区间,这个范围的合理选择将大大增加找到最佳轨迹的机会。

(2)飞行状态2

边界约束:φ(tf)=φ(t0)+π/2,xi(tf)=xi(t0),yi(tf)=yi(t0),飞行轨迹近似为圆形,也称为游荡轨迹。φ的范围是[-π,π].

(3)飞行状态3

边界约束:φ(tf)=φ(t0)+π,xi(tf)和yi(tf)无约束。 这个飞行轨迹是基本轨迹的一个变体,初始和最终航向方向正好相反,即航向方向的变化180°,φ的范围是[-π/2,3π/2]。

三种飞行状态的最佳飞行轨迹如图2所示。图2中的轨迹与三种飞行状态的边界条件完全一致。从图2中可以看出,高度的范围很大,这是因为有一种假设,即风场的强度随着海拔高度的增加而增加。由于无人机体积较大,它需要飞得更高才能获得更多的风能。无人机的初始速度也很大,因此前方飞行距离也比通常情况下大。

图2 三种飞行状态的最佳飞行轨迹

在上述三种飞行状态中,系统的状态变量为[xi,yi,h,Vt,γ,φ],控制变量为[C1,μ]。[xi,yi,h,Vt]、[γ,φ]、[C1,μ]的变化曲线和系统的能量变化分别如图3(a)～5(a)、3(b)～5(b)、3(c)～5(c)和3(d)～5(d)所示。通过比较图3至图5中的曲线,可以看到真空速Vt在三种状态下变化趋势相同, 变化范围为40～130 m·s-1。爬升角/下滑角γ在-50°到50°之间,γ的变化趋势相同。φ的范围为边界约束,分别是π/2,2π和π。升力系数CL在飞行状态1中达到最大值,大约为1.3,而CL在飞行状态3中小于0.9,飞行状态2中的CL大小在其他两种状态之间。μ的变化趋势和CL相似,μ大部分时间保持正值,这意味着无人机以正方向飞行,这与图2一致,图2中的三个轨迹都沿着x轴的正方向增加。此外μ被限幅在±75°之间且μ在三种飞行状态下都达到极限值并保持一段时间。 三种飞行状态的能量变化趋势大致相同,重力势能先增大后减小,动能先减小后增大,总能量趋于恒定。

图3 飞行状态1的变量曲线

图4 飞行状态2的变量曲线

图5 飞行状态3的变量曲线

3.2 固定风场下的最小功率滑翔

考虑具有恒定推力的动力飞行,公式(23)可改写为

(49)

图6 和最小推力T的关系

图7 三种取值情况下的最优轨迹

图8 恒定推力下的各变量曲线

4 结论

本文提出了一种无人机动态翱翔飞行轨迹优化方法,给出了获得无人机动态翱翔运动方程的具体步骤;选择了梯度风场模型和惯性坐标系,使这些方程与大多数文献的表达式一致;采用Radau伪谱方法对飞行轨迹进行优化。由于所选无人机比普通无人机大,因此在仿真中需要正确选择参数。仿真结果表明,该方法可以成功地生成无人机动态飞行的最优轨迹。此外,通过对无动力飞行和有动力飞行的比较,发现最佳发动机推力和风场强度近似成反比。