CN104573236A - Continuous system simulation verification-based partial finite element method - Google Patents

Abstract

The invention discloses a continuous system simulation verification-based partial finite element method. The method aims at the defect that a dynamic optimization system is solved by a current combination method, and the continuous system simulation verification-based partial finite element method is provided by only guaranteeing that the constraint on a finite element configuration point is met but not guaranteeing that the constraint on a non-configuration point is met. According to the method, the constraint is met on the finite element non-configuration point by using a continuous system simulation verification method; rearrangement is performed on an original finite element sequence by using the partial finite element method, so a finite element configuration scheme can guarantee the solving accuracy on the premise that the number of finite elements is increased as less as possible. Solving is performed by using the method aiming at a free flight conflict resolution problem in a three-dimensional space; the solving result shows that the discretization precision and the solving precision are met better by using the method compared with the original combination method.

Description

A kind of subdivision Finite Element Method based on continuous system simulation checking

Technical field

The invention belongs to the track optimizing field with end conswtraint and profile constraints, relate to a kind of subdivision Finite Element Method based on continuous system simulation checking.

Background technology

Under current air traffic control pattern, the contradiction between the air traffic day by day increased and limited spatial domain resource.Therefore, free flight becomes the inevitable development pattern of following air traffic.Under free flight environment, user can from the main separation line of flight, flying speed and flying height, need not by the air route of a series of guidance station, course line fetter; And can according to the optimization aim of oneself as saving fuel oil, shorten boat time, dodge hazardous weather or spatial domain and hinder, determine flight path, thus reach making full use of spatial domain resource.Free flight will bring the advantage of spatial domain high power capacity, but simultaneously under free flight condition, to increase considerably in the same time, with the aircraft density occurred in spatial domain, add the uncertainty in free course line, thus cause the possibility of flight collision significantly to promote.For ensureing being perfectly safe of flight, when predict will clash time, need to design the Optimum Operation strategy avoiding flight collision, make aircraft according to the optimal trajectory flight of setting, thus break away from the flight collision that may occur.Therefore, how devise optimum operation strategy realizes Solving Flight Conflicts is affect the gordian technique and condition precedent that can free flight realize.

The present invention carries out dynamic optimization research to conflict Resolution problem between the multi-aircraft in the stipulated time, dynamics mass-point model is all adopted to describe for each airplane, and the initial position of multi-aircraft and final position are fixed, all need the constraint strictly meeting Air Traffic System separation standard between every two airplanes, proposition optimization aim is that aircraft consumption gross energy is minimum.

At present, simultaneous method is the main stream approach solving optimization problems.Control variable and state variable are carried out sliding-model control by simultaneous method in whole optimization time domain simultaneously, original optimization problems is changed into large-scale nonlinear programming (NLP) problem, the nonlinear programming problem after then utilizing nonlinear planning solution device to solve discretize.This temporal partitioning is multiple time period by finite element by simultaneous method, and the connection coordination by Law based on orthogonal collocation on finite element meets the collocation point of interpolation point as finite element of orthogonal polynomial.On each collocation point of finite element, differential algebraic equations is strict satisfied.The solution of subordination principle (DAE) system combines with optimization problem by simultaneous method, only needs to solve a DAE system of equations at optimum solution place, so just avoids the calculating that pilot process is a large amount of.In addition, simultaneous method can process path constraint problem very easily, has good stability simultaneously.Nonlinear planning solution device of the present invention is the IPOPT developed by CMU.

When adopting simultaneous method to solve dynamic optimization proposition, it is more intensive that finite element divides, the precision of the solution tried to achieve is higher, the scale of nonlinear programming problem is larger, the solution procedure of nonlinear planning solution device is more consuming time more difficult, therefore the number only by increasing finite element ensures that the precision of separating is infeasible at some time, especially for the extensive optimization problems of complexity.The more important thing is, when adopting simultaneous method to solve dynamic optimization proposition at present, can only ensure that limited first wife puts a little upper satisfied constraint, not ensure non-collocation point meets constraint.This defect is reflected in free flight problem, namely show as this proposition and be optimized the flight path solving and successfully obtain, seem between aircraft and all meet constraint, actually on the non-collocation point of some finite element, because thus the constraint not meeting Air Traffic System separation standard there occurs conflict, this just causes the solving result of simultaneous method to allow people relieved completely.The present invention is directed to this defect that simultaneous method solves dynamic optimization proposition, namely simultaneous method is merely able to ensure to meet constraint on finite element collocation point, do not ensure to meet constraint on non-collocation point, propose a kind of subdivision Finite Element Method based on continuous system simulation checking, continuous system simulation verification method makes the non-collocation point of finite element meet constraint, subdivision Finite Element Method reconfigures original limit metasequence, can obtain reduce as far as possible add finite element number prerequisite under ensure the finite element allocation plan of solving precision, solving result shows to adopt to solve optimal problem in this way and well meets discretize precision and solving precision.

Summary of the invention

The present invention is directed to the defect of full simultaneous discretize dynamic optimization at present, namely simultaneous method is merely able to ensure to meet constraint on finite element collocation point, do not ensure to meet constraint on non-collocation point, propose a kind of subdivision Finite Element Method based on continuous system simulation checking, can reduce as far as possible add finite element number prerequisite under ensure the finite element allocation plan of solving precision.

The conflict Resolution in Free Flight problem that the present invention is directed under three dimensions is studied, state Solving Flight Conflicts problem as optimum control proposition, simultaneous method is adopted to solve it, and with the collision error of continuous system simulation checking for foundation, subdivision Finite Element Method is adopted to reconfigure finite element, to improve the precision solved, after trying to achieve control variable curve, continuous simulation is carried out by integrator, and compare with the numerical result before finite element reconfigures, the validity of verification algorithm.

Method of the present invention comprises the following steps:

Step (1). to flight collision time domain t to be measured _fbe divided into Nfe section, Nfe≤5, obtain Nfe finite element, then the finite element sequence after waiting point time domain is each finite element length is

Step (2). adopt orthogonal collocation on finite element method to be configured the finite element sequence in step (1), have 3 Radau collocation points in each finite element, the point between every two collocation points is non-collocation point, specifically:

The general type of dynamic optimization proposition is:

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(t\right),p)\≤0\\ t\∈[t_0,t_f]\end{array}\right.---\left(1\right)$

Wherein φ objective function, differential variable, algebraic variable, control variable, model parameter, n _z, n _y, n _u, n _pthe number of differential variable, algebraic variable, control variable, model parameter respectively, z ₀it is the original state of differential variable.F represents the differential equation, g _erepresentation algebra equation, in subordination principle, g is summed up in the point that in the boundary constraint of differential variable, algebraic variable and control variable _iin.

Orthogonal collocation on finite element method is by the orthonormal polynomial approximation control variable in finite element and state variable, and definition finite element number is NFe, so t ₀< t ₁< ... < t _nfe=t _f, h _i=t _i-t _i-1.Each finite element h _ion by Lagrange interpolation polynomial, differential variable, control variable and algebraic variable are approached:

$\left\{\begin{array}{c}z\left(t\right)={\mathrm{\&Sigma;}}_{j=0}^K{l}_j\left(\mathrm{\&tau;}\right)z_{i,j}\\ u\left(t\right)={\mathrm{\&Sigma;}}_{j=1}^K\stackrel{\&OverBar;}{l_J}\left(\mathrm{\&tau;}\right)u_{i,j}\\ y\left(t\right)={\mathrm{\&Sigma;}}_{j=1}^K\stackrel{\&OverBar;}{l_J}\left(\mathrm{\&tau;}\right)y_{i,j}\\ t=t_{i-1}+h_i\mathrm{\&tau;}\end{array}\right.t\&Element;[t_{i-1},t_i],\mathrm{\&tau;}\&Element;\left[\mathrm{0,1}\right]---\left(2\right)$

Wherein, K is the order of interpolation, l _j(τ) and represent the Lagrange interpolation polynomial of differential variable and control variable, algebraic variable respectively, can represent by following form:

$\left\{\begin{array}{c}{l}_j\left(\mathrm{\&tau;}\right)={\mathrm{\&Pi;}}_{k=0,k\&NotEqual;j}^K\frac{\mathrm{\&tau;}-{\mathrm{\&tau;}}_k}{{\mathrm{\&tau;}}_j-{\mathrm{\&tau;}}_k}\\ \stackrel{\&OverBar;}{l_J}\left(\mathrm{\&tau;}\right)={\mathrm{\&Pi;}}_{k=1,k\&NotEqual;j}^K\frac{\mathrm{\&tau;}-{\mathrm{\&tau;}}_k}{{\mathrm{\&tau;}}_j-{\mathrm{\&tau;}}_k}\end{array}\right.---\left(3\right)$

Lagrange polynomial interpolation has following character:

$l_j\left({\mathrm{\&tau;}}_k\right)=\left\{\begin{array}{c}1,k=j\\ 0,k\&NotEqual;j\end{array}\right.,\stackrel{\&OverBar;}{l_J}\left({\mathrm{\&tau;}}_k\right)=\left\{\begin{array}{c}1,k=j\\ 0,k\&NotEqual;j\end{array}\right.---\left(4\right)$

Namely the value of each variable on collocation point just in time equals its coefficient, so

z(t _i，j)＝z _i，j，y(t _i，j)＝y _i，j，u(t _i，j)＝u _i，j(5)

Because differential variable needs the continuity of hold mode, so need the continuity being ensured differential variable by Connection equations on finite element end points, algebraic variable and control variable then can be discontinuous.In addition, initial and terminal condition must also be added:

$z_{i+1}={\mathrm{\&Sigma;}}_{j=0}^K{l}_j\left(1\right)z_{i,j},z_{\mathrm{1,0}}=z\left(t_0\right),z_f={\mathrm{\&Sigma;}}_{j=0}^K{l}_j\left(1\right)z_{\mathrm{NE},j}---\left(6\right)$

Formula (2) and (4) are substituted into (1) and in conjunction with condition of contact equation (6), after can obtaining discretize, NLP problem form is as follows:

${\min }_{z_{i,j},y_{i,j},u_{i,j}}\mathrm{\&phi;}\left(z\left(t_f\right)\right)$

$s.t.\left\{\begin{array}{c}{\mathrm{\&Sigma;}}_{k=0}^K{l}_k^{\&prime;}\left({\mathrm{\&tau;}}_j\right)z_{i,k}-h_i\&CenterDot;f(z_{i,j},y_{i,j},u_{i,j},p)=0,i=1,...,\mathrm{Nfe};j=1,...,K\\ g_I(z_{i,j},y_{i,j},u_{i,j},p)\&le;0,i=1,...,\mathrm{Nfe};j=1,...,K\\ g_E(z_{i,j},y_{i,j},u_{i,j},p)=0,i=1,...,\mathrm{Nfe};j=1,...,K\\ z_{i+\mathrm{1,0}}={\mathrm{\&Sigma;}}_{j=0}^K{l}_j\left(1\right)z_{i,j},i=1,...,\mathrm{Nfe}-1\\ z_{\mathrm{1,0}}=z\left(t_0\right),z_f={\mathrm{\&Sigma;}}_{j=0}^K{l}_j\left(1\right)z_{\mathrm{Nfe},j}\end{array}\right.---\left(7\right)$

Step (3). the NLP problem after the discretize obtain step 2 and formula (7) solve, and obtain Optimal Flight Route scheme:

The present invention investigates the conflict Resolution in Free Flight problem of N frame aircraft.Range of movement due to research object aircraft far exceedes itself physical dimension, therefore by abstract for aircraft be particle.For each frame aircraft, following kinetic model is adopted to be described,

$\left\{\begin{array}{c}\frac{{\mathrm{dx}}_t\left(t\right)}{\mathrm{dt}}=v_{x_t}\left(t\right)\\ \frac{{\mathrm{dy}}_t\left(t\right)}{\mathrm{dt}}=v_{y_t}\left(t\right)\\ \frac{{\mathrm{dz}}_t\left(t\right)}{\mathrm{dt}}=v_{z_t}\left(t\right)\end{array}\right.---\left(8\right)$

i＝1....N

Wherein, (x _i(t), y _i(t), z _i(t)) be the position coordinates of the i-th airplane in t, unit is m;

be the speed of the i-th frame aircraft in t, unit is m/s, and it is the control variable handling aircraft.

According to the flight safety boundary condition of ATC (air traffic control) standard, any two horizontal range Rs of frame aircraft both when sustained height _ijbe not less than R (R=5nmi) or both vertical range H _ijbe not less than H (H=1000ft), following expression formula of extracting can be adopted to be described:

Wherein, R _ijbe the horizontal range between the i-th frame aircraft and jth frame aircraft, unit is m;

H _ijvertical range between i-th airplane and jth frame aircraft, unit is m.

Assuming that need at t ₀to t _fcomplete the process of conflict Resolution during this period of time, original state and the final state of each frame aircraft be,

$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)$

Wherein, be the initial position co-ordinates of the i-th frame aircraft, unit is m;

be the final position coordinate of the i-th frame aircraft, unit is m.

T ₀for starting the moment of carrying out conflict Resolution, unit is s, t _ffor terminating the moment of conflict Resolution, unit is s, and t _f-t ₀=15min

In general, the good and bad performance index weighing each frame aircraft running orbit of conflict Resolution process mainly contain time, fuel consumption or energy ezpenditure etc.Case definition objective function of the present invention is,

$J={\mathrm{\&Sigma;}}_{i=1}^N{v}_i\frac{1}{2}{\&Integral;}_{t=t_0}^{t_f}{({\left(\frac{{\mathrm{dx}}_t}{\mathrm{dt}}\right)}^2+{\left(\frac{{\mathrm{dy}}_t}{\mathrm{dt}}\right)}^2+{\mathrm{\&eta;}}^2{\left(\frac{{\mathrm{dz}}_t}{\mathrm{dt}}\right)}^2)}^2\mathrm{dt}---\left(11\right)$

Wherein, J is that objective function is namely using the energy ezpenditure in Solving Flight Conflicts process as performance index;

η is penalty factor, because speed change in vertical direction can cause the discomfort of airborne personnel, therefore adopts η to punish the velocity variations on z direction;

V _ifor the weight coefficient of each airplane.

In sum, with the speed component v of N airplane _xi, v _yi, v _zifor control variable, the Solving Flight Conflicts problem of N airplane is expressed as the optimum control proposition of following form,

$\left\{\begin{array}{c}\min J\\ s.t.\mathrm{Equ}\left(1\right)~\mathrm{Equ}\left(3\right)\end{array}\right.---\left(12\right)$

For formula (2) in the relational expression (R that extracts _ij>=R) ∨ (H _ij>=H), the people such as Raghunathan are by introducing auxiliary variable λ _ij1and λ _ij2be translated into following form,

$\left\{\begin{array}{c}{\mathrm{\&lambda;}}_{\mathrm{ij}1}(R_{\mathrm{ij}}-R)+{\mathrm{\&lambda;}}_{\mathrm{ij}2}(H_{\mathrm{ij}}-H)\&GreaterEqual;0\\ {\mathrm{\&lambda;}}_{\mathrm{ij}1}+{\mathrm{\&lambda;}}_{\mathrm{ij}2}=1\\ {\mathrm{\&lambda;}}_{\mathrm{ij}1}\&GreaterEqual;0,{\mathrm{\&lambda;}}_{\mathrm{ij}2}\&GreaterEqual;0\end{array}\right.---\left(13\right)$

Introduce about λ in the objective function of optimum control proposition simultaneously _ij1and λ _ij2quadratic term,

$A_{\mathrm{ij}}={\&Integral;}_{t=t_0}^{t_f}{\left({\mathrm{\&lambda;}}_{\mathrm{ij}1}\right)}^2+{\left({\mathrm{\&lambda;}}_{\mathrm{ij}2}\right)}^2\mathrm{dt}---\left(14\right)$

Therefore optimum control proposition becomes following form,

$\left\{\begin{array}{c}\min J^{\&prime;}=J+\frac{1}{2}{\mathrm{\&Sigma;}}_{i=1}^N{\mathrm{\&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)$

By the air speed component v of the objective function J ' correspondence in formula (8) _xi, v _yi, v _ziconstitute Optimal Flight Route scheme, also carry out continuous system simulation checking for step (3) simultaneously and provide control variable.

Step (4). continuous system simulation checking is carried out to the non-collocation point of finite element that step 2 obtains, obtains the flight path scheme of non-collocation point:

4.1: the control variable v obtained according to step (3) Optimization Solution _xi, v _yi, v _zi, adopt Matlab this control variable effect to be put in continuous system and carry out simulating, verifying, can state variable x be obtained _i, y _i, z _ivalue on the non-collocation point of finite element, i.e. the flight path scheme of the non-collocation point of finite element;

The present invention adopts the ode45 () function in matlab to carry out simulating, verifying to continuous system.

4.2: for the flight path scheme of the non-collocation point of finite element that step (4.1) obtains, on each non-collocation point, whether every two airplanes are clashed (see formula 16) according to the secure border condition judgment between aircraft in whole flight course; If clash, then find the non-collocation point that this conflict is corresponding, perform step (5); If no conflict occurred, then proposition Optimization Solution success, namely obtains Optimal Flight Route scheme.

From the secure border condition between aircraft, the horizontal range of any two airplanes both when sustained height is not less than to be thought when R or both vertical ranges are not less than H and does not collide.Therefore, in flight course, two airplanes clash and meet formula (9) i.e. horizontal range be less than R and vertical range be less than H,

Wherein, x _i(t _noc), y _i(t _noc), z _i(t _noc) represent that the i-th airplane is at t _nocthe position coordinates in moment (the non-collocation point of finite element)

Step (5). according to the non-collocation point t clashed found in step (4) _noc, adopt subdivision Finite Element Method at this non-collocation point, obtain a new finite element sequence, return and perform step (2).

Subdivision finite element is that the finite element corresponding to this non-collocation point carries out subdivision, then obtains new finite element sequence.Suppose there be Nfe finite element at present, finite element sequence is verify that the non-collocation point clashed obtained is t according to continuous system simulation _noc, corresponding M the finite element of this point.Therefore, the new finite element sequence obtained through subdivision Finite Element Method is

Beneficial effect of the present invention is as follows:

The present invention is directed to the defect of full simultaneous discretize dynamic optimization at present, namely simultaneous method is merely able to ensure to meet constraint on finite element collocation point, do not ensure to meet constraint on non-collocation point, propose a kind of subdivision Finite Element Method based on continuous system simulation checking, distich legislation improves, continuous system simulation verification method makes the non-collocation point of finite element meet constraint, subdivision Finite Element Method reconfigures original limit metasequence, can reduce as far as possible add finite element number prerequisite under ensure the finite element allocation plan of solving precision.

The present invention is directed to the conflict Resolution in Free Flight problem under three dimensions, state Solving Flight Conflicts problem as optimum control proposition, employing continuous system simulation verification method obtains the conflict moment (non-collocation point) in flight course, subdivision Finite Element Method is adopted to reconfigure finite element sequence, solving result shows to adopt to solve optimal problem in this way and well meets discretize precision and solving precision, avoid defect when former simultaneous method solves optimum control proposition, demonstrate the feasibility of the method.

Accompanying drawing explanation

Fig. 1 is the process flow diagram of the method for the invention operational process.

Embodiment

Below in conjunction with drawings and Examples, the present invention is described in further detail.

The course of work of the method for the invention as shown in Figure 1,

Step (1). to flight collision time domain t to be measured _fbe divided into Nfe section, Nfe≤5, obtain Nfe finite element, then the finite element sequence after waiting point time domain is each finite element length is

Step (2). adopt orthogonal collocation on finite element method to be configured the finite element sequence in step (1), have 3 Radau collocation points in each finite element, the point between every two collocation points is non-collocation point:

The general type of dynamic optimization proposition is:

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(t\right),p)\&le;0\\ t\&Element;[t_0,t_f]\end{array}\right.---\left(18\right)$

Wherein φ objective function, differential variable, algebraic variable, control variable, model parameter, n _z, n _y, n _u, n _pthe number of differential variable, algebraic variable, control variable, model parameter respectively, z ₀it is the original state of differential variable.F represents the differential equation, g _erepresentation algebra equation, in subordination principle, g is summed up in the point that in the boundary constraint of differential variable, algebraic variable and control variable _iin.

Orthogonal collocation on finite element method is by the orthonormal polynomial approximation control variable in finite element and state variable.Definition finite element number is NFe, so t ₀< t ₁< ... < t _nfe=t _f, h _i=t _i-t _i-1.Each finite element h _ion by Lagrange interpolation polynomial, differential variable, control variable and algebraic variable are approached:

$\left\{\begin{array}{c}z\left(t\right)={\mathrm{\&Sigma;}}_{j=0}^K{l}_j\left(\mathrm{\&tau;}\right)z_{i,j}\\ u\left(t\right)={\mathrm{\&Sigma;}}_{j=1}^K\stackrel{\&OverBar;}{l_J}\left(\mathrm{\&tau;}\right)u_{i,j}\\ y\left(t\right)={\mathrm{\&Sigma;}}_{j=1}^K\stackrel{\&OverBar;}{l_J}\left(\mathrm{\&tau;}\right)y_{i,j}\\ t=t_{i-1}+h_i\mathrm{\&tau;}\end{array}\right.t\&Element;[t_{i-1},t_i],\mathrm{\&tau;}\&Element;\left[\mathrm{0,1}\right]---\left(19\right)$

Wherein, K is the order of interpolation, l _j(τ) and represent the Lagrange interpolation polynomial of differential variable and control variable, algebraic variable respectively, can represent by following form:

$\left\{\begin{array}{c}{l}_j\left(\mathrm{\&tau;}\right)={\mathrm{\&Pi;}}_{k=0,k\&NotEqual;j}^K\frac{\mathrm{\&tau;}-{\mathrm{\&tau;}}_k}{{\mathrm{\&tau;}}_j-{\mathrm{\&tau;}}_k}\\ \stackrel{\&OverBar;}{l_J}\left(\mathrm{\&tau;}\right)={\mathrm{\&Pi;}}_{k=1,k\&NotEqual;j}^K\frac{\mathrm{\&tau;}-{\mathrm{\&tau;}}_k}{{\mathrm{\&tau;}}_j-{\mathrm{\&tau;}}_k}\end{array}\right.---\left(20\right)$

Lagrange polynomial interpolation has following character:

$l_j\left({\mathrm{\&tau;}}_k\right)=\left\{\begin{array}{c}1,k=j\\ 0,k\&NotEqual;j\end{array}\right.,\stackrel{\&OverBar;}{l_J}\left({\mathrm{\&tau;}}_k\right)=\left\{\begin{array}{c}1,k=j\\ 0,k\&NotEqual;j\end{array}\right.---\left(21\right)$

Namely the value of each variable on collocation point just in time equals its coefficient, so

z(t _i，j)＝z _i，j，y(t _i，j)＝y _i，j，u(t _i，j)＝u _i，j(22)

Because differential variable needs the continuity of hold mode, so need the continuity being ensured differential variable by Connection equations on finite element end points, algebraic variable and control variable then can be discontinuous.In addition, initial and terminal condition must also be added:

$z_{i+1}={\mathrm{\&Sigma;}}_{j=0}^K{l}_j\left(1\right)z_{i,j},z_{\mathrm{1,0}}=z\left(t_0\right),z_f={\mathrm{\&Sigma;}}_{j=0}^K{l}_j\left(1\right)z_{\mathrm{NE},j}---\left(23\right)$

Formula (19) and (21) are substituted into (18) and in conjunction with condition of contact equation (23), after can obtaining discretize, NLP problem form is as follows:

${\min }_{z_{i,j},y_{i,j},u_{i,j}}\mathrm{\&phi;}\left(z\left(t_f\right)\right)$

$s.t.\left\{\begin{array}{c}{\mathrm{\&Sigma;}}_{k=0}^K{l}_k^{\&prime;}\left({\mathrm{\&tau;}}_j\right)z_{i,k}-h_i\&CenterDot;f(z_{i,j},y_{i,j},u_{i,j},p)=0,i=1,...,\mathrm{Nfe};j=1,...,K\\ g_I(z_{i,j},y_{i,j},u_{i,j},p)\&le;0,i=1,...,\mathrm{Nfe};j=1,...,K\\ g_E(z_{i,j},y_{i,j},u_{i,j},p)=0,i=1,...,\mathrm{Nfe};j=1,...,K\\ z_{i+\mathrm{1,0}}={\mathrm{\&Sigma;}}_{j=0}^K{l}_j\left(1\right)z_{i,j},i=1,...,\mathrm{Nfe}-1\\ z_{\mathrm{1,0}}=z\left(t_0\right),z_f={\mathrm{\&Sigma;}}_{j=0}^K{l}_j\left(1\right)z_{\mathrm{Nfe},j}\end{array}\right.---\left(24\right)$

Step (3). the NLP problem after the discretize obtain step 2 solves, and obtains Optimal Flight Route scheme:

The present invention investigates the conflict Resolution in Free Flight problem of N frame aircraft.Range of movement due to research object aircraft far exceedes itself physical dimension, therefore by abstract for aircraft be particle.For each frame aircraft, following kinetic model is adopted to be described,

$\left\{\begin{array}{c}\frac{{\mathrm{dx}}_t\left(t\right)}{\mathrm{dt}}=v_{x_t}\left(t\right)\\ \frac{{\mathrm{dy}}_t\left(t\right)}{\mathrm{dt}}=v_{y_t}\left(t\right)\\ \frac{{\mathrm{dz}}_t\left(t\right)}{\mathrm{dt}}=v_{z_t}\left(t\right)\end{array}\right.---\left(25\right)$

i＝1....N

Wherein, (x _i(t), y _i(t), z _i(t)) be the position coordinates of the i-th airplane in t, unit is m;

be the speed of the i-th frame aircraft in t, unit is m/s, and it is the control variable handling aircraft.

According to the flight safety boundary condition of ATC (air traffic control) standard, any two horizontal range Rs of frame aircraft both when sustained height _ijbe not less than R (R=5nmi) or both vertical range H _ijbe not less than H (H=1000ft), following expression formula of extracting can be adopted to be described:

Wherein, R _ijbe the horizontal range between the i-th frame aircraft and jth frame aircraft, unit is m;

H _ijbe the vertical range between the i-th airplane and jth frame aircraft, unit is m.

Assuming that need at t ₀to t _fcomplete the process of conflict Resolution during this period of time, original state and the final state of each frame aircraft be,

$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(27\right)$

Wherein, be the initial position co-ordinates of the i-th frame aircraft, unit is m;

be the final position coordinate of the i-th frame aircraft, unit is m.

T ₀for starting the moment of carrying out conflict Resolution, unit is s, t _ffor terminating the moment of conflict Resolution, unit is s, and t _f-t ₀=15min

In general, the good and bad performance index weighing each frame aircraft running orbit of conflict Resolution process mainly contain time, fuel consumption or energy ezpenditure etc.Case definition objective function of the present invention is,

$J={\mathrm{\&Sigma;}}_{i=1}^N{v}_i\frac{1}{2}{\&Integral;}_{t=t_0}^{t_f}{({\left(\frac{{\mathrm{dx}}_t}{\mathrm{dt}}\right)}^2+{\left(\frac{{\mathrm{dy}}_t}{\mathrm{dt}}\right)}^2+{\mathrm{\&eta;}}^2{\left(\frac{{\mathrm{dz}}_t}{\mathrm{dt}}\right)}^2)}^2\mathrm{dt}---\left(28\right)$

Wherein, J is that objective function is namely using the energy ezpenditure in Solving Flight Conflicts process as performance index;

η is penalty factor, because speed change in vertical direction can cause the discomfort of airborne personnel, therefore adopts η to punish the velocity variations on z direction;

V _ifor the weight coefficient of each airplane.

In sum, with the speed component v of N airplane _xi, v _yi, v _zifor control variable, the Solving Flight Conflicts problem of N airplane is expressed as the optimum control proposition of following form,

$\left\{\begin{array}{c}\min J\\ s.t.\mathrm{Equ}\left(25\right)~\mathrm{Equ}\left(27\right)\end{array}\right.---\left(29\right)$

For formula (2) in the relational expression (R that extracts _ij>=R) ∨ (H _ij>=H), the people such as Raghunathan are by introducing auxiliary variable λ _ij1and λ _ij2be translated into following form,

$\left\{\begin{array}{c}{\mathrm{\&lambda;}}_{\mathrm{ij}1}(R_{\mathrm{ij}}-R)+{\mathrm{\&lambda;}}_{\mathrm{ij}2}(H_{\mathrm{ij}}-H)\&GreaterEqual;0\\ {\mathrm{\&lambda;}}_{\mathrm{ij}1}+{\mathrm{\&lambda;}}_{\mathrm{ij}2}=1\\ {\mathrm{\&lambda;}}_{\mathrm{ij}1}\&GreaterEqual;0,{\mathrm{\&lambda;}}_{\mathrm{ij}2}\&GreaterEqual;0\end{array}\right.---\left(30\right)$

Introduce about λ in the objective function of optimum control proposition simultaneously _ij1and λ _ij2quadratic term,

$A_{\mathrm{ij}}={\&Integral;}_{t=t_0}^{t_f}{\left({\mathrm{\&lambda;}}_{\mathrm{ij}1}\right)}^2+{\left({\mathrm{\&lambda;}}_{\mathrm{ij}2}\right)}^2\mathrm{dt}---\left(31\right)$

Therefore optimum control proposition becomes following form,

$\left\{\begin{array}{c}\min J^{\&prime;}=J+\frac{1}{2}{\mathrm{\&Sigma;}}_{i=1}^N{\mathrm{\&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(32\right)$

By the air speed component v of the objective function J ' correspondence in formula (8) _xi, v _yi, v _ziconstitute Optimal Flight Route scheme, also carry out continuous system simulation checking for step (3) simultaneously and provide control variable.

Step (4). continuous system simulation checking is carried out to the non-collocation point of finite element that step 2 obtains, obtains the flight path scheme of non-collocation point:

4.1: the control variable v obtained according to step (3) Optimization Solution _xi, v _yi, v _zi, adopt Matlab this control variable effect to be put in continuous system and carry out simulating, verifying, can state variable x be obtained _i, y _i, z _ivalue on the non-collocation point of finite element, i.e. the flight path scheme of the non-collocation point of finite element;

The present invention adopts the ode45 () function in matlab to carry out simulating, verifying to continuous system.

4.2: for the flight path scheme of the non-collocation point of finite element that step (4.1) obtains, on each non-collocation point, whether every two airplanes are clashed (see formula 33) according to the secure border condition judgment between aircraft in whole flight course; If clash, then find the non-collocation point that this conflict is corresponding, perform step (5); If no conflict occurred, then proposition Optimization Solution success, namely obtains Optimal Flight Route scheme.

From the secure border condition between aircraft, the horizontal range of any two airplanes both when sustained height is not less than to be thought when R or both vertical ranges are not less than H and does not collide.Therefore, in flight course, two airplanes clash and meet formula (9) i.e. horizontal range be less than R and vertical range be less than H,

Wherein, x _i(t _noc), y _i(t _noc), z _i(t _noc) represent that the i-th airplane is at t _nocthe position coordinates in moment (the non-collocation point of finite element)

Step (5). according to the non-collocation point t clashed found in step (4) _noc, adopt subdivision Finite Element Method at this non-collocation point, obtain a new finite element sequence, return and perform step (2).

Subdivision finite element is that the finite element corresponding to this non-collocation point carries out subdivision, then obtains new finite element sequence.Suppose there be Nfe finite element at present, finite element sequence is verify that the non-collocation point clashed obtained is t according to continuous system simulation _noc, corresponding M the finite element of this point.Therefore, the new finite element sequence obtained through subdivision Finite Element Method is

Claims (4)

1., based on a subdivision Finite Element Method for continuous system simulation checking, it is characterized in that the method comprises the following steps:

Step (1). to flight collision time domain t to be measured _fbe divided into Nfe section, Nfe≤5, obtain Nfe finite element, then the finite element sequence after waiting point time domain is each finite element length is

Step (2). adopt orthogonal collocation on finite element method to be configured the finite element sequence in step (1), have 3 Radau collocation points in each finite element, the point between every two collocation points is non-collocation point, specifically:

By the conflict Resolution in Free Flight problem statement order already issued topic under three dimensions, the general type of dynamic optimization proposition is:

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\&Element;[t_0,t_f]\end{array}\right.---\left(1\right);$

Wherein φ objective function, differential variable, algebraic variable, control variable, model parameter, n _z, n _y, n _u, n _pthe number of differential variable, algebraic variable, control variable, model parameter respectively, z ₀it is the original state of differential variable.F represents the differential equation, g _brepresentation algebra equation, in subordination principle, g is summed up in the point that in the boundary constraint of differential variable, algebraic variable and control variable _iin;

Orthogonal collocation on finite element method is by the orthonormal polynomial approximation control variable in finite element and state variable, and definition finite element number is Nfe, so t ₀< t ₁< ... < t _nfe=t _f, h _i=t _i-t _i-1; Each finite element h _ion by Lagrange interpolation polynomial, differential variable, control variable and algebraic variable are approached:

$\left\{\begin{array}{c}z\left(t\right)={\mathrm{\&Sigma;}}_{j=0}^K{l}_j\left(\mathrm{\&tau;}\right)z_{i,j}\\ u\left(t\right)={\mathrm{\&Sigma;}}_{j=1}^K\stackrel{\&OverBar;}{l_J}\left(\mathrm{\&tau;}\right)u_{i,j}\\ y\left(t\right)={\mathrm{\&Sigma;}}_{j=1}^K\stackrel{\&OverBar;}{l_J}\left(\mathrm{\&tau;}\right)y_{i,j}\\ t=t_{i-1}+h_i\mathrm{\&tau;}\end{array}\right.t\&Element;[t_{i-1},t_i],\mathrm{\&tau;}\&Element;\left[\mathrm{0,1}\right]---\left(2\right);$

Wherein, K is the order of interpolation, l _j(τ) and represent differential variable and control variable respectively, the Lagrange interpolation polynomial of algebraic variable, can represent by following form:

$\left\{\begin{array}{c}{l}_j\left(\mathrm{\&tau;}\right)={\mathrm{\&Pi;}}_{k=0,k\&NotEqual;j}^K\frac{\mathrm{\&tau;}-{\mathrm{\&tau;}}_k}{{\mathrm{\&tau;}}_j-{\mathrm{\&tau;}}_k}\\ \stackrel{\&OverBar;}{l_J}\left(\mathrm{\&tau;}\right)={\mathrm{\&Pi;}}_{k=1,k\&NotEqual;j}^K\frac{\mathrm{\&tau;}-{\mathrm{\&tau;}}_k}{{\mathrm{\&tau;}}_j-{\mathrm{\&tau;}}_k}\end{array}\right.---\left(3\right);$

Lagrange polynomial interpolation has following character:

$l_j\left({\mathrm{\&tau;}}_k\right)=\left\{\begin{array}{c}1,k=j\\ 0,k\&NotEqual;j\end{array}\right.,\stackrel{\&OverBar;}{l_J}\left({\mathrm{\&tau;}}_k\right)=\left\{\begin{array}{c}1,k=j\\ 0,k\&NotEqual;\end{array}\right.---\left(4\right);$

Namely the value of each variable on collocation point just in time equals its coefficient, so

z(t _i，j)＝z _i，j，y(t _i，j)＝y _i，j，u(t _i，j)＝u _i，j(5)；

Because differential variable needs the continuity of hold mode, so need the continuity being ensured differential variable by Connection equations on finite element end points, algebraic variable and control variable then can be discontinuous; In addition, initial and terminal condition must also be added:

$z_{i+1}={\mathrm{\&Sigma;}}_{j=0}^K{l}_j\left(1\right)z_{i,j},z_{\mathrm{1,0}}=z\left(t_0\right),z_f={\mathrm{\&Sigma;}}_{j=0}^K{l}_j\left(1\right)z_{\mathrm{NE},j}---\left(6\right);$

Formula (2) and (4) are substituted into (1) and in conjunction with condition of contact equation (6), after can obtaining discretize, NLP problem form is as follows:

Step (3). the NLP problem after the discretize obtain step 2 solves, and obtains Optimal Flight Route scheme:

Following kinetic model is adopted to be described each frame aircraft:

$\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

Wherein, (x _i(t), y _i(t), z _i(t)) be the position coordinates of the i-th airplane in t, unit is m;

be the speed of the i-th frame aircraft in t, unit is m/s, and it is the control variable handling aircraft;

According to the flight safety boundary condition of ATC standard, any two horizontal range Rs of frame aircraft both when sustained height _ijbe not less than R (R=5nmi) or both vertical range H _ijbe not less than H (H=1000ft), following expression formula of extracting can be adopted to be described:

Wherein, R _ijbe the horizontal range between the i-th frame aircraft and jth frame aircraft, unit is m;

H _ijbe the vertical range between the i-th airplane and jth frame aircraft, unit is m;

Assuming that need at t ₀to t _fcomplete the process of conflict Resolution during this period of time, original state and the final state of each frame aircraft are:

$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);$

Wherein, be the initial position co-ordinates of the i-th frame aircraft, unit is m;

be the final position coordinate of the i-th frame aircraft, unit is m;

T ₀for starting the moment of carrying out conflict Resolution, unit is s;

T _ffor terminating the moment of conflict Resolution, unit is s, and t _f-t ₀=15min;

This method is using the energy ezpenditure in Solving Flight Conflicts process as performance index, and namely objective function J is shown in formula (11):

$j={\mathrm{\&Sigma;}}_{i=1}^N{v}_i\frac{1}{2}{\&Integral;}_{t=t_0}^{t_f}{({\left(\frac{{\mathrm{dx}}_i}{\mathrm{dt}}\right)}^2+{\left(\frac{{\mathrm{dy}}_i}{\mathrm{dt}}\right)}^2+{\mathrm{\&eta;}}^2{\left(\frac{{\mathrm{dz}}_i}{\mathrm{dt}}\right)}^2)}^2\mathrm{dt}---\left(11\right);$

Wherein η is penalty factor, because speed change in vertical direction can cause the discomfort of airborne personnel, therefore adopts η to punish the velocity variations on z direction; v _ifor the weight coefficient of each airplane;

In sum, with the speed component v of N airplane _xi, v _yi, v _zifor control variable, the Solving Flight Conflicts problem of N airplane is expressed as the optimum control proposition of following form:

$\left\{\begin{array}{c}\min J\\ s.t.\mathrm{Equ}\left(1\right)~\mathrm{Equ}\left(3\right)\end{array}\right.---\left(12\right);$

For formula (2) in relational expression of extracting will by introducing auxiliary variable λ _ij1and λ _ij2be translated into following form:

$\left\{\begin{array}{c}{\mathrm{\&lambda;}}_{\mathrm{ij}1}(R_{\mathrm{ij}}-R)+{\mathrm{\&lambda;}}_{\mathrm{ij}2}(H_{\mathrm{ij}}-H)\&GreaterEqual;0\\ {\mathrm{\&lambda;}}_{\mathrm{ij}1}+{\mathrm{\&lambda;}}_{\mathrm{ij}2}=1\\ {\mathrm{\&lambda;}}_{\mathrm{ij}1}\&GreaterEqual;0,{\mathrm{\&lambda;}}_{\mathrm{ij}2}\&GreaterEqual;0\end{array}\right.---\left(13\right);$

Introduce about λ in the objective function of optimum control proposition simultaneously _ij1and λ _ij2quadratic term:

$A_{\mathrm{ij}}={\&Integral;}_{t=t_0}^{t_f}{\left({\mathrm{\&lambda;}}_{\mathrm{ij}1}\right)}^2+{\left({\mathrm{\&lambda;}}_{\mathrm{ij}2}\right)}^2\mathrm{dt}---\left(14\right);$

Therefore optimum control proposition becomes following form:

$\left\{\begin{array}{c}\min J^{\&prime;}=J+\frac{1}{2}{\mathrm{\&Sigma;}}_{i=1}^N{\mathrm{\&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);$

By the air speed component v of the objective function I ' correspondence in formula (8) _xi, v _yi, v _ziconstitute Optimal Flight Route scheme, also carry out continuous system simulation checking for step (3) simultaneously and provide control variable;

Step (4). continuous system simulation checking is carried out to the non-collocation point of finite element that step (2) obtains, obtains the flight path scheme of non-collocation point:

Step (5). according to the non-collocation point t clashed found in step (4) _noc, adopt subdivision Finite Element Method at this non-collocation point, obtain a new finite element sequence, return and perform step (2).

2. as claimed in claim 1 a kind of based on continuous system simulation checking subdivision Finite Element Method, it is characterized in that step (4) concrete operations are as follows:

4.1: the control variable v obtained according to step (3) Optimization Solution _xi, v _yi, v _zi, adopt Matlab this control variable effect to be put in continuous system and carry out simulating, verifying, can state variable x be obtained _i, y _i, z _ivalue on the non-collocation point of finite element, i.e. the flight path scheme of the non-collocation point of finite element;

4.2: for the flight path scheme of the non-collocation point of finite element that step (4.1) obtains, on each non-collocation point, whether every two airplanes are clashed according to the secure border condition judgment between aircraft in whole flight course, sees formula 16; If clash, then find the non-collocation point that this conflict is corresponding, perform step (5); If no conflict occurred, then proposition Optimization Solution success, namely obtains Optimal Flight Route scheme;

In described flight course, two airplanes clash, then demand fulfillment formula (9) i.e. horizontal range be less than R and vertical range be less than H:

Wherein, x _i(t _noc), y _i(t _noc), z _i(t _noc) represent that the i-th airplane is at t _nocthe position coordinates in moment, t _nocbe the non-collocation point of finite element.

3. a kind of subdivision Finite Element Method based on continuous system simulation checking as claimed in claim 2, is characterized in that step (4.1) adopts the ode45 () function in matlab to carry out simulating, verifying to continuous system.

4. as claimed in claim 2 a kind of based on continuous system simulation checking subdivision Finite Element Method, it is characterized in that the subdivision finite element described in step (5) is that the finite element corresponding to this non-collocation point carries out subdivision, then new finite element sequence is obtained, concrete operations suppose there be Nfe finite element at present, and finite element sequence is verify that the non-collocation point clashed obtained is t according to continuous system simulation _noc, corresponding M the finite element of this point.Therefore, the new finite element sequence obtained through subdivision Finite Element Method is: