一种基于连续系统仿真验证的剖分有限元方法

摘要

本发明公开一种基于连续系统仿真验证的剖分有限元方法。该方法针对目前联立法求解动态优化命题的缺陷，即只保证在有限元配置点上满足约束，不保证在非配置点上满足约束，提出了一种基于连续系统仿真验证的剖分有限元方法。该方法采用连续系统仿真验证方法使得有限元非配置点满足约束，采用剖分有限元方法对原有限元序列进行重新配置，能够得到在尽量少增加有限元个数的前提下保证求解精度的有限元配置方案。本发明针对三维空间下的自由飞行冲突解脱问题采用该方法进行求解，求解结果表明，与原联立法相比，采用这种方法更好的满足了离散化精度和求解精度。

主权项

一种基于连续系统仿真验证的剖分有限元方法，其特征在于该方法包括以下步骤：步骤(1).对待测飞行冲突时域t_f等分为Nfe段，Nfe≤5，得到Nfe个有限元，则经过等分时域后的有限元序列为每个有限元长度为步骤(2).采用有限元正交配置方法对步骤(1)中的有限元序列进行配置，每个有限元中有3个Radau配置点，每两个配置点之间的点均为非配置点，具体是：将三维空间下的自由飞行冲突解脱问题表述成命题，动态优化命题的一般形式为：minφ(z(t_f))$\left\{\begin{array}{c}\frac{\mathrm{dz}\left(t\right)}{\mathrm{dt}}=f(z\left(t\right),y\left(t\right),u\left(t\right),p),z\left(0\right)=z_0\\ g_E(z\left(t\right),y\left(t\right),u\left(t\right),p)=0\\ g_I(z\left(t\right),y\left(t\right),u\left(y\right),p)\le 0\\ t\in [t_0,t_f]\end{array}\right.---\left(1\right);$其中是φ目标函数，是微分变量，是代数变量，是控制变量，是模型参数，n_z、n_y、n_u、n_p分别是微分变量、代数变量、控制变量、模型参数的个数，z₀是微分变量的初始状态。f表示微分方程，g_B表示代数方程，微分代数方程组中微分变量、代数变量和控制变量的边界约束归结到g_I中；有限元正交配置方法通过有限元上的正交多项式逼近控制变量和状态变量，定义有限元个数为Nfe，那么t₀＜t₁＜…＜t_Nfe＝t_f，h_i＝t_i‑t_i‑1；每个有限元h_i上通过Lagrange插值多项式对微分变量、控制变量和代数变量进行逼近：$\left\{\begin{array}{c}z\left(t\right)={\Sigma }_{j=0}^K{l}_j\left(\tau \right)z_{i,j}\\ u\left(t\right)={\Sigma }_{j=1}^K\overline{l_J}\left(\tau \right)u_{i,j}\\ y\left(t\right)={\Sigma }_{j=1}^K\overline{l_J}\left(\tau \right)y_{i,j}\\ t=t_{i-1}+h_i\tau \end{array}\right.t\in [t_{i-1},t_i],\tau \in \left[\mathrm{0,1}\right]---\left(2\right);$其中，K为插值的阶次，l_j(τ)和分别表示微分变量和控制变量,代数变量的Lagrange插值多项式，可以用如下形式表示：$\left\{\begin{array}{c}{l}_j\left(\tau \right)={\Pi }_{k=0,k\ne j}^K\frac{\tau -{\tau }_k}{{\tau }_j-{\tau }_k}\\ \overline{l_J}\left(\tau \right)={\Pi }_{k=1,k\ne j}^K\frac{\tau -{\tau }_k}{{\tau }_j-{\tau }_k}\end{array}\right.---\left(3\right);$Lagrange多项式插值具有如下性质：$l_j\left({\tau }_k\right)=\left\{\begin{array}{c}1,k=j\\ 0,k\ne j\end{array}\right.,\overline{l_J}\left({\tau }_k\right)=\left\{\begin{array}{c}1,k=j\\ 0,k\ne \end{array}\right.---\left(4\right);$即各个变量在配置点上的值正好等于其系数，那么z(t_i，j)＝z_i，j，y(t_i，j)＝y_i，j，u(t_i，j)＝u_i，j    (5)；由于微分变量需要保持状态的连续性，所以在有限元端点上需要通过连接方程来保证微分变量的连续性，而代数变量和控制变量则可以不连续；此外，还必须加上初始和终端条件：$z_{i+1}={\Sigma }_{j=0}^K{l}_j\left(1\right)z_{i,j},z_{\mathrm{1,0}}=z\left(t_0\right),z_f={\Sigma }_{j=0}^K{l}_j\left(1\right)z_{\mathrm{NE},j}---\left(6\right);$将公式(2)和(4)代入(1)并结合连接条件方程(6)，可以得到离散化后NLP问题形式如下：步骤(3).对步骤2得到的离散化后的NLP问题进行求解，得到最优飞行轨迹方案：采用如下的动力学模型进行描述每一架飞行器：$\left\{\begin{array}{c}\frac{{\mathrm{dx}}_i\left(t\right)}{\mathrm{dt}}=v_{x_i}\left(t\right)\\ \frac{{\mathrm{dy}}_i\left(t\right)}{\mathrm{dt}}=v_{y_i}\left(t\right)\\ \frac{{\mathrm{dz}}_i\left(t\right)}{\mathrm{dt}}=v_{z_i}\left(t\right)\end{array}\right.---\left(8\right);$i＝1....N其中，(x_i(t)，y_i(t)，z_i(t))为第i架飞机在t时刻的位置坐标，单位为m；为第i架飞行器在t时刻的速度，单位为m/s，且是操纵飞行器的控制变量；根据ATC标准的飞行安全边界条件可知，任意两架飞行器在同一高度时两者的水平距离R_ij不小于R(R＝5nmi)或两者的垂直距离H_ij不小于H(H＝1000ft)，可以采用如下的析取表达式进行描述：其中，R_ij为第i架飞行器和第j架飞行器之间的水平距离，单位为m；H_ij为第i架飞机和第j架飞行器之间的垂直距离，单位为m；假定需在t₀到t_f这段时间内完成冲突解脱的过程，各架飞行器的初始状态和终止状态为：$x_i\left(0\right)=x_{i,t_0},y_i\left(0\right)=y_{i,t_0},z_i\left(0\right)=z_{i,t_0}$$x_i\left(t_f\right)=x_{i,t_f},y_i\left(t_f\right)=y_{i,t_f},z_i\left(t_f\right)=z_{i,t_f}---\left(10\right);$其中，为第i架飞行器的初始位置坐标，单位为m；为第i架飞行器的终止位置坐标，单位为m；t₀为开始进行冲突解脱的时刻，单位为s；t_f为结束冲突解脱的时刻，单位为s,且t_f‑t₀＝15min；本方法将飞行冲突解脱过程中的能量消耗作为性能指标，即目标函数J见公式(11)：$j={\Sigma }_{i=1}^N{v}_i\frac{1}{2}{\int }_{t=t_0}^{t_f}{({\left(\frac{{\mathrm{dx}}_i}{\mathrm{dt}}\right)}^2+{\left(\frac{{\mathrm{dy}}_i}{\mathrm{dt}}\right)}^2+{\eta }^2{\left(\frac{{\mathrm{dz}}_i}{\mathrm{dt}}\right)}^2)}^2\mathrm{dt}---\left(11\right);$其中η为惩罚因子，由于在垂直方向上的速度改变会引起机载人员的不适，因此采用η对z方向上的速度变化进行惩罚；v_i为各架飞机的权重系数；综上所述，以N架飞机的速度分量v_xi，v_yi，v_zi为控制变量，将N架飞机的飞行冲突解脱问题表述为如下形式的最优控制命题：$\left\{\begin{array}{c}\min J\\ s.t.\mathrm{Equ}\left(1\right)~\mathrm{Equ}\left(3\right)\end{array}\right.---\left(12\right);$针对式⑵中的析取关系式将通过引入辅助变量λ_ij1和λ_ij2将其转化为如下形式：$\left\{\begin{array}{c}{\lambda }_{\mathrm{ij}1}(R_{\mathrm{ij}}-R)+{\lambda }_{\mathrm{ij}2}(H_{\mathrm{ij}}-H)\ge 0\\ {\lambda }_{\mathrm{ij}1}+{\lambda }_{\mathrm{ij}2}=1\\ {\lambda }_{\mathrm{ij}1}\ge 0,{\lambda }_{\mathrm{ij}2}\ge 0\end{array}\right.---\left(13\right);$同时在最优控制命题的目标函数中引入关于λ_ij1和λ_ij2的二次项：$A_{\mathrm{ij}}={\int }_{t=t_0}^{t_f}{\left({\lambda }_{\mathrm{ij}1}\right)}^2+{\left({\lambda }_{\mathrm{ij}2}\right)}^2\mathrm{dt}---\left(14\right);$故最优控制命题变为如下形式：$\left\{\begin{array}{c}\min J^{\prime }=J+\frac{1}{2}{\Sigma }_{i=1}^N{\Sigma }_{j=i+1}^N{A}_{\mathrm{ij}}\\ s.t.\mathrm{Equ}\left(8\right),\mathrm{Equ}\left(10\right),\mathrm{Equ}\left(13\right)\end{array}\right.---\left(15\right);$由公式(8)中的目标函数I′对应的飞机速度分量v_xi，v_yi，v_zi构成了最优飞行轨迹方案，同时也为步骤(3)进行连续系统仿真验证提供了控制变量；步骤(4).对步骤(2)得到的有限元非配置点进行连续系统仿真验证，得到非配置点的飞行轨迹方案：步骤(5).根据步骤(4)中找到的发生冲突的非配置点t_noc，在该非配置点采用剖分有限元方法，得到一个新的有限元序列，返回执行步骤(2)。