Collision-free Multiple Unmanned Combat Aerial Vehicles Cooperative Trajectory Planning for Time-critical Missions using Differential Flatness Approach

This paper investigates the cooperative trajectory planning for multiple unmanned combat aerial vehicles in performing autonomous cooperative air-to-ground target attack missions. Firstly, the collision-free cooperative trajectory planning problem for time-critical missions is formulated as a cooperative trajectory optimal control problem (CTP-OCP), which is based on an approximate allowable attack region model, several constraints model, and a multi-criteria objective function. Next, a planning algorithm based on the differential flatness, B-spline curves and nonlinear programming is designed to solve the CTP-OCP. In particular, the notion of the virtual time is introduced to deal with the temporal constraints. Finally, the proposed approach is validated by two typical scenarios and the simulation results show the feasibility and effectiveness of the proposed planning approach.


Keywords:    Cooperative trajectory planning,   unmanned combat aerial vehicle,   differential flatness,  B-spline curve,  allowable attack region 

Nowadays, it is an active research area to perform autonomous cooperative air-to-ground target attack (CA/GTA) missions by multiple unmanned combat aerial vehicles (multi-UCAV)1 . However, it is quite difficult to coordinate the strike operation of multi-UCAV, especially for time-critical missions that require precise timing and sequencing of tasks and operations. In order to accomplish the mission successfully, UCAVs need to generate detailed cooperative execution plans to lead themselves well to penetrate through enemy threat areas, avoid the collisions with obstacles or each other, fly into the allowable attack regions (AARs) simultaneously or in sequence, and then perform weapon delivery operations. The cooperative trajectory planning for CA/GTA is vital to achieve the mission goals. It is really one of the key challenging technologies, due to its high dimensionality, severe equality and inequality constraints involved, and the requirement of spatial-temporal cooperation of multi-UCAV, and has recently received extensive attentions2 .

To date, various algorithms3 - 7 have been developed to solve this cooperative planning problem, including several collision-avoidance techniques and time adjustment strategies. McLain3 , et al. used coordination variables and coordination functions based strategies to achieve cooperative timing among teams of vehicles, by coordinating the velocity and path length of each vehicle. Kaminer4 , et al. proposed a solution to the problem of coordinated control of multiple unmanned aerial vehicles (multi-UAV) to ensure collision-free maneuvers under strict spatial and temporal constraints. Bollino5 , et al. addressed the optimal path planning of UAVs in obstacle-rich environments and proposed the collision-free path planning for multi-UAV using optimal control techniques and pseudospectral methods. Lian6 introduced a differential flatness based approach to optimally formulate the cooperative path planning for multi-agent dynamical systems considering spatial and temporal constraints, and parameterized the curves by B-spline representations. Kuwata7 , et al. presented a cooperative distributed robust trajectory optimization approach, using RH-MILP with independent dynamics but coupled objectives and hard constraints. The above investigations have given several valuable strategies in the cooperative planning, but they failed to tackle the point-to-region cooperative trajectory planning for CA/GTA missions under consideration directly, which needs to integrate both the spatial and temporal constraints on the level of the trajectory planning. To address the trajectory planning, a novel cooperative trajectory planning algorithm for multi-UCAV in performing the CA/GTA missions is presented

Given a set of stationary ground targets in a terrain region, the mission objective is to plan several cooperative optimal or suboptimal, dynamically feasible flight trajectories from the initial points (IPs) into the AARs, such that multi-UCAV can effectively attack the targets with minimum time, maximum survivability and target destruction effectiveness. Therefore, many factors need to be considered synthetically in the cooperative trajectory planning. In this section, the cooperative trajectory planning problem is formulated after the modelling of AAR, constraints and the objective function.

2.1  Allowable Attack Region Model


For the point-to-region trajectory planning problem, the AARs of targets are defined as the areas where UCAVs can effectively perform weapon delivery operations. So, in order to plan the accurate and optimal attack trajectories for weapon delivery, the AARs and the delivery parameters need to be obtained. The AAR of the ith target TARi, denoted as R(TARi), is such a set of all feasible release states that TARi can be effectively attacked whenever the aircraft is in that states. R(TARi) can be formulated as an abstract 6-dimensional space8


R(TA R i )={ x,y,h,v,γ,ψ } 6 . MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaGaamOuaiaacI cacaWGubGaamyqaiaadkfadaWgaaWcbaGaamyAaaqabaGccaGGPaGa eyypa0ZaaiWaaeaacaWG4bGaaiilaiaadMhacaGGSaGaamiAaiaacY cacaWG2bGaaiilaiabeo7aNjaacYcacqaHipqEaiaawUhacaGL9baa cqGHckcZcqWIDesOdaahaaWcbeqaaiaaiAdaaaGccaqGUaaaaa@507F@       (1)

Obviously, R(TARi) is high-dimensional nonlinear space, which will make the problem solving much difficult. By presetting an appropriate weapon release speed, vr and the flight-path angle, γ r MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabmqaamaabaabaaGcbaGaeq4SdC2aaS baaSqaaiaabkhaaeqaaaaa@38A0@ based on estimating the weapon impact effects and destruction requirements to the target, and predetermining release heading ψ r MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaGafqiYdKNbau aadaWgaaWcbaGaaeOCaaqabaaaaa@3A9D@ ,it can be reduced to a 3-dimensional space


R(TA R i )LA R sp (TA R i ,[ v r , γ r ]) ={ (x,y,h, v r , γ r , ψ r )| (x,y,h)LAR(TA R i , [ v r , γ r , ψ r ]) ψ r =azimuth of  TA R i relative to (x,y,h) }. MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGceaqabeaacaWGsb GaaiikaiaadsfacaWGbbGaamOuamaaBaaaleaacaWGPbaabeaakiaa cMcacqGHshI3caWGmbGaamyqaiaadkfadaWgaaWcbaGaae4Caiaabc haaeqaaOGaaiikaiaadsfacaWGbbGaamOuamaaBaaaleaacaWGPbaa beaakiaacYcacaGGBbGaamODamaaBaaaleaacaqGYbaabeaakiaacY cacqaHZoWzdaWgaaWcbaGaaeOCaaqabaGccaGGDbGaaiykaaqaaiab g2da9maacmaabaGaaiikaiaadIhacaGGSaGaamyEaiaacYcacaWGOb GaaiilaiaadAhadaWgaaWcbaGaaeOCaaqabaGccaGGSaGaeq4SdC2a aSbaaSqaaiaabkhaaeqaaOGaaiilaiqbeI8a5zaafaWaaSbaaSqaai aabkhaaeqaaOGaaiykamaaeeaaeaqabeaacaGGOaGaamiEaiaacYca caWG5bGaaiilaiaadIgacaGGPaGaeyicI4SaamitaiaadgeacaWGsb GaaiikaiaadsfacaWGbbGaamOuamaaBaaaleaacaWGPbaabeaakiaa cYcacaqGGaGaai4waiaadAhadaWgaaWcbaGaaeOCaaqabaGccaGGSa Gaeq4SdC2aaSbaaSqaaiaabkhaaeqaaOGaaiilaiqbeI8a5zaafaWa aSbaaSqaaiaabkhaaeqaaOGaaiyxaiaacMcaaeaacuaHipqEgaqbam aaBaaaleaacaqGYbaabeaakiabg2da9iaabggacaqG6bGaaeyAaiaa b2gacaqG1bGaaeiDaiaabIgacaqGGaGaae4BaiaabAgacaqGGaGaae iiaiaadsfacaWGbbGaamOuamaaBaaaleaacaWGPbaabeaakiaamcca caqGYbGaaeyzaiaabYgacaqGHbGaaeiDaiaabMgacaqG2bGaaeyzai aabccacaqG0bGaae4BaiaabccacaGGOaGaamiEaiaacYcacaWG5bGa aiilaiaadIgacaGGPaaaaiaawEa7aaGaay5Eaiaaw2haaiaac6caaa aa@A482@       (2)

2.2  Constraint Model


2.2.1 Maneuverability Constraints of UCAVs

The maneuverability constraints impact every phase during target attack missions, so they should be met for an executable plan. The constraints can be denoted as



v min v(t) v max ,      γ min γ(t) γ max ,       ψ min ψ(t) ψ max , μ min μ(t) μ max ,   n x,min n x (t) n x,max ,  n h,min n h (t) n h,max , MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGceaqabeaacaWG2b WaaSbaaSqaaiGac2gacaGGPbGaaiOBaaqabaGccqGHKjYOcaWG2bGa aiikaiaadshacaGGPaGaeyizImQaamODamaaBaaaleaaciGGTbGaai yyaiaacIhaaeqaaOGaaiilaiaabccacaqGGaGaaeiiaiaabccacaqG GaGaeq4SdC2aaSbaaSqaaiGac2gacaGGPbGaaiOBaaqabaGccqGHKj YOcqaHZoWzcaGGOaGaamiDaiaacMcacqGHKjYOcqaHZoWzdaWgaaWc baGaciyBaiaacggacaGG4baabeaakiaacYcacaqGGaGaaeiiaiaabc cacaqGGaGaaeiiaiaabccacqaHipqEdaWgaaWcbaGaciyBaiaacMga caGGUbaabeaakiabgsMiJkabeI8a5jaacIcacaWG0bGaaiykaiabgs MiJkabeI8a5naaBaaaleaaciGGTbGaaiyyaiaacIhaaeqaaOGaaiil aaqaaiabeY7aTnaaBaaaleaaciGGTbGaaiyAaiaac6gaaeqaaOGaey izImQaeqiVd0MaaiikaiaadshacaGGPaGaeyizImQaeqiVd02aaSba aSqaaiGac2gacaGGHbGaaiiEaaqabaGccaGGSaGaaeiiaiaabccaca WGUbWaaSbaaSqaaiaabIhacaGGSaGaciyBaiaacMgacaGGUbaabeaa kiabgsMiJkaad6gadaWgaaWcbaGaaeiEaaqabaGccaGGOaGaamiDai aacMcacqGHKjYOcaWGUbWaaSbaaSqaaiaabIhacaGGSaGaciyBaiaa cggacaGG4baabeaakiaacYcacaqGGaGaamOBamaaBaaaleaacaqGOb GaaiilaiGac2gacaGGPbGaaiOBaaqabaGccqGHKjYOcaWGUbWaaSba aSqaaiaabIgaaeqaaOGaaiikaiaadshacaGGPaGaeyizImQaamOBam aaBaaaleaacaqGObGaaeilaiaab2gacaqGHbGaaeiEaaqabaGccaGG Saaaaaa@AC83@       (3)

where v is the aircraft velocity; γ MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaGaeq4SdCgaaa@394A@ denotes the flight-path angle; denotes the heading angle ψ MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaGaeqiYdKhaaa@3971@ ; denotes the roll angle μ MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaGaeqiVd0gaaa@3959@ ; nx and nh are the longitudinal and normal components of the load factor. and () min MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaGaaiikaiabgw SixlaacMcadaWgaaWcbaGaciyBaiaacMgacaGGUbaabeaaaaa@3E43@ and () max MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaGaaiikaiabgw SixlaacMcadaWgaaWcbaGaciyBaiaacggacaGG4baabeaaaaa@3E45@ represent the minimum and maximum boundary values of the value () MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaGaaiikaiabgw SixlaacMcaaaa@3B45@ .


2.2.2 Battlefield Environment Constraints

(a) Threat avoidance: without loss of generality, the threat impact range can be assumed as hemispheres approximately, thus threat avoidance constraint can be defined as


( x(t) x i o ,y(t) y i o ,h(t) h i o ) 2 R i o , (i=1,2,, N o ), MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaWaauWaaeaada qadaqaaiaadIhacaGGOaGaamiDaiaacMcacqGHsislcaWG4bWaa0ba aSqaaiaadMgaaeaacaqGVbaaaOGaaiilaiaadMhacaGGOaGaamiDai aacMcacqGHsislcaWG5bWaa0baaSqaaiaadMgaaeaacaqGVbaaaOGa aiilaiaadIgacaGGOaGaamiDaiaacMcacqGHsislcaWGObWaa0baaS qaaiaadMgaaeaacaqGVbaaaaGccaGLOaGaayzkaaaacaGLjWUaayPc SdWaaSbaaSqaaiaaikdaaeqaaOGaeyyzImRaamOuamaaDaaaleaaca WGPbaabaGaae4BaaaakiaacYcacaqGGaGaaeikaiaadMgacqGH9aqp caaIXaGaaiilaiaaikdacaGGSaGaeS47IWKaaiilaiaad6eadaWgaa WcbaGaae4BaaqabaGccaGGPaGaaiilaaaa@65E9@       (4)

where the norm 2 MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipy0lf9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabeqaamaabaabaaGcbaWaauWaaKaaGf aacqGHflY1aiaawMa7caGLkWoakmaaBaaajeaybaGaaGOmaaqabaaa aa@3CA3@ denotes the Euclidean distance between two points and ( x i o , y i o , h i o ) MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaGaaiikaiaadI hadaqhaaWcbaGaamyAaaqaaiaab+gaaaGccaGGSaGaamyEamaaDaaa leaacaWGPbaabaGaae4BaaaakiaacYcacaWGObWaa0baaSqaaiaadM gaaeaacaqGVbaaaOGaaiykaaaa@4389@ and R i o MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaGaamOuamaaDa aaleaacaWGPbaabaGaae4Baaaaaaa@3A87@ denote the origin coordinates and the detection radius of the ith threat, respectively.
(b) No-fly zones (NFZs): herein, an infinite-length cylinder is used to describe the NFZs:

( x(t) x i w ,y(t) y i w ) 2 R i w , (i=1,2,, N w ), MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaWaauWaaeaada qadaqaaiaadIhacaGGOaGaamiDaiaacMcacqGHsislcaWG4bWaa0ba aSqaaiaadMgaaeaacaqG3baaaOGaaiilaiaadMhacaGGOaGaamiDai aacMcacqGHsislcaWG5bWaa0baaSqaaiaadMgaaeaacaqG3baaaaGc caGLOaGaayzkaaaacaGLjWUaayPcSdWaaSbaaSqaaiaaikdaaeqaaO GaeyyzImRaamOuamaaDaaaleaacaWGPbaabaGaae4DaaaakiaacYca caqGGaGaaeikaiaadMgacqGH9aqpcaaIXaGaaiilaiaaikdacaGGSa GaeS47IWKaaiilaiaad6eadaWgaaWcbaGaae4DaaqabaGccaGGPaGa aiilaaaa@5E29@       (5)

where ( x i w , y i w ) MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaGaaiikaiaadI hadaqhaaWcbaGaamyAaaqaaiaabEhaaaGccaGGSaGaamyEamaaDaaa leaacaWGPbaabaGaae4DaaaakiaacMcaaaa@3FE4@ and R i w MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaGaamOuamaaDa aaleaacaWGPbaabaGaae4Daaaaaaa@3A8F@ denote the origin coordinates and the radius of the ith NFZ.

2.2.3 Terminal Constraints

In the point-to-region trajectory planning, weapon delivery point (WDPt) as the trajectory terminal of CA/GTA needs to be in the AAR, i.e. satisfying the terminal constraints. The formula can be denoted as WDPt i R(TA R i ) MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipy0lf9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabeqaamaabaabaaGcbaGaae4vaiaabs eacaqGqbGaaeiDamaaBaaaleaacaWGPbaabeaakiabgIGiolaadkfa caqGOaGaamivaiaadgeacaWGsbWaaSbaaSqaaiaadMgaaeqaaOGaae ykaaaa@4153@ , that is


| x a x( t f ) |Δx| y a y( t f ) |Δy, | h a h( t f ) |Δh, MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaWaaqWaaeaaca WG4bWaaSbaaSqaaiaabggaaeqaaOGaeyOeI0IaamiEaiaacIcacaWG 0bWaaSbaaSqaaGqaaiaa=zgaaeqaaOGaaiykaaGaay5bSlaawIa7ai abgsMiJkabfs5aejaadIhacaqGSaGaaeiiamaaemaabaGaamyEamaa BaaaleaacaqGHbaabeaakiabgkHiTiaadMhacaGGOaGaamiDamaaBa aaleaacaWFMbaabeaakiaacMcaaiaawEa7caGLiWoacqGHKjYOcqqH uoarcaWG5bGaaiilaiaabccadaabdaqaaiaadIgadaWgaaWcbaGaae yyaaqabaGccqGHsislcaWGObGaaiikaiaadshadaWgaaWcbaGaa8Nz aaqabaGccaGGPaaacaGLhWUaayjcSdGaeyizImQaeuiLdqKaamiAai aacYcaaaa@66D1@       (6)

where tf is the terminal time, (xa, ya, ha) and (x(tf), y(tf), h(tf)) are the coordinates of AAR’s center point and WDPt, (Δx, Δy, Δh) MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaqcaaSaaiikaK aaGkabfs5aejaadIhajaaWcaGGSaGaaeiiaKaaGkabfs5aejaadMha caGGSaGaaeiiaiabfs5aejaadIgajaaWcaGGPaaaaa@4669@ are the thresholds of errors.

2.2.4 Cooperative Constraints

During the mission, multi-UCAV should maintain a safe distance to avoid collision with each other, i.e. satisfying spatial constraint. The model can be denoted as


ρ j ( t i ) ρ k ( t i ) 2 max( d safe j ,  d safe k ), ( j,k=1,2,, N v ,jk), MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaWaauWaaeaacq aHbpGCdaahaaWcbeqaaiaadQgaaaGccaGGOaGaamiDamaaBaaaleaa caWGPbaabeaakiaacMcacqGHsislcqaHbpGCdaahaaWcbeqaaiaadU gaaaGccaGGOaGaamiDamaaBaaaleaacaWGPbaabeaakiaacMcaaiaa wMa7caGLkWoadaWgaaWcbaGaaGOmaaqabaGccqGHLjYSciGGTbGaai yyaiaacIhacaGGOaGaamizamaaDaaaleaacaqGZbGaaeyyaiaabAga caqGLbaabaGaamOAaaaakiaacYcacaqGGaGaamizamaaDaaaleaaca qGZbGaaeyyaiaabAgacaqGLbaabaGaam4AaaaakiaacMcacaGGSaGa aeiiaiaabIcacaqGGaGaamOAaiaacYcacaWGRbGaeyypa0JaaGymai aacYcacaaIYaGaaiilaiabl+UimjaacYcacaWGobWaaSbaaSqaaiaa bAhaaeqaaOGaaiilaiaadQgacqGHGjsUcaWGRbGaaiykaiaacYcaaa a@6F7C@       (7)

where ρ j ( t i )={ x j ( t i ), y j ( t i ), h j ( t i )} MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaGaeqyWdi3aaW baaSqabeaacaWGQbaaaOGaaiikaiaadshadaWgaaWcbaGaamyAaaqa baGccaGGPaGaeyypa0Jaai4EaiaadIhadaahaaWcbeqaaiaadQgaaa GccaGGOaGaamiDamaaBaaaleaacaWGPbaabeaakiaacMcacaGGSaGa amyEamaaCaaaleqabaGaamOAaaaakiaacIcacaWG0bWaaSbaaSqaai aadMgaaeqaaOGaaiykaiaacYcacaWGObWaaWbaaSqabeaacaWGQbaa aOGaaiikaiaadshadaWgaaWcbaGaamyAaaqabaGccaGGPaGaaiyFaa aa@5320@ is the spatial position of UCAVj at the time ti,   d safe j MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaGaaeiiaiaads gadaqhaaWcbaGaae4CaiaabggacaqGMbGaaeyzaaqaaiaadQgaaaaa aa@3DF5@ is the minimum safety radius of UCAVj, and Nv is the total number of UCAVs.

The temporal constraints considered here include simultaneous arrival constraint and tight-sequencing constraint, which can be described as

t f j + Δ jk t f k 0, (j,k=1,2,, N v ,jk), MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaGaamiDamaaDa aaleaacaqGMbaabaGaamOAaaaakiabgUcaRiabfs5aenaaBaaaleaa caWGQbGaam4AaaqabaGccqGHsislcaWG0bWaa0baaSqaaiaabAgaae aacaWGRbaaaOGaeyizImQaaGimaiaacYcacaqGGaGaaeikaiaadQga caGGSaGaam4Aaiabg2da9iaaigdacaGGSaGaaGOmaiaacYcacqWIVl ctcaGGSaGaamOtamaaBaaaleaacaqG2baabeaakiaacYcacaWGQbGa eyiyIKRaam4AaiaacMcacaGGSaaaaa@5829@       (8)

where t f j MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipy0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabeqaamaabaabaaGcbaGaamiDamaaDa aaleaacaqGMbaabaGaamOAaaaaaaa@38D5@ is the terminal time of UCAVj, Δ jk MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaGaeuiLdq0aaS baaSqaaiaadQgacaWGRbaabeaaaaa@3B13@ is the arrival interval between UCAVs.

2.3   Objective Function


The objective function of each UCAVj can be defined by the weighted sum of the three separate running cost terms with appropriate weighting factors w t j MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaGaam4DamaaDa aaleaacaqG0baabaGaamOAaaaaaaa@3AB1@ w r j MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaGaam4DamaaDa aaleaacaqGYbaabaGaamOAaaaaaaa@3AAF@ and w d j MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaGaam4DamaaDa aaleaacaqGKbaabiqaamtbcaWGQbaaaaaa@3ADC@ , where w t j + w r j + w d j =1 MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaGaam4DamaaDa aaleaacaqG0baabaGaamOAaaaakiabgUcaRiaadEhadaqhaaWcbaGa aeOCaaqaaiaadQgaaaGccqGHRaWkcaWG3bWaa0baaSqaaiaabsgaae aacaWGQbaaaOGaeyypa0JaaGymaaaa@4460@ . And the objective function of the entire team can defined as:


J=min j=1 N v J j =min j=1 N v ( w t j J t j + w r j J prd j + w d j J dest j ) , MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaGaamOsaiabg2 da9iGac2gacaGGPbGaaiOBamaaqahabaGaamOsamaaCaaaleqabaGa amOAaaaaaeaacaWGQbGaeyypa0JaaGymaaqaaiaad6eadaWgaaadba GaamODaaqabaaaniabggHiLdGccqGH9aqpciGGTbGaaiyAaiaac6ga daaeWbqaaiaacIcacaWG3bWaa0baaSqaaiaabshaaeaacaWGQbaaaO GaamOsamaaDaaaleaacaqG0baabaGaamOAaaaakiabgUcaRiaadEha daqhaaWcbaGaaeOCaaqaaiaadQgaaaGccaWGkbWaa0baaSqaaiaabc hacaqGYbGaaeizaaqaaiaadQgaaaGccqGHRaWkcaWG3bWaa0baaSqa aiaabsgaaeaacaWGQbaaaOGaamOsamaaDaaaleaacaqGKbGaaeyzai aabohacaqG0baabaGaamOAaaaakiaacMcaaSqaaiaadQgacqGH9aqp caaIXaaabaGaamOtamaaBaaameaacaWG2baabeaaa0GaeyyeIuoaki aacYcaaaa@6A56@          (9)

and three separate running cost terms can be respectively defined as

J t j = 1 Γ j t 0 j t f j dt = 1 Γ j ( t f j t 0 j ), MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaGaamOsamaaDa aaleaacaqG0baabaGaamOAaaaakiabg2da9maalaaabaGaaGymaaqa aiabfo5ahnaaCaaaleqabaGaamOAaaaaaaGcdaWdXaqaaiaadsgaca WG0baaleaacaWG0bWaa0baaWqaaiaaicdaaeaacaWGQbaaaaWcbaGa amiDamaaDaaameaacaqGMbaabaGaamOAaaaaa0Gaey4kIipakiabg2 da9maalaaabaGaaGymaaqaaiabfo5ahnaaCaaaleqabaGaamOAaaaa aaGccaGGOaGaamiDamaaDaaaleaacaqGMbaabaGaamOAaaaakiabgk HiTiaadshadaqhaaWcbaGaaGimaaqaaiaadQgaaaGccaGGPaGaaiil aaaa@5612@          (10)

J prd j 1 t f j t 0 j t 0 j t f j PRD j (t)dt = 1 t f j t 0 j t 0 j t f j [1 r=1 n r (1 P d j (t,r)) ]dt , MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaGaamOsamaaDa aaleaacaqGWbGaaeOCaiaabsgaaeaacaWGQbaaaOGaaeypaiaabcca daWcaaqaaiaaigdaaeaacaWG0bWaa0baaSqaaiaabAgaaeaacaWGQb aaaOGaeyOeI0IaamiDamaaDaaaleaacaaIWaaabaGaamOAaaaaaaGc daWdXaqaaiaabcfacaqGsbGaaeiramaaCaaaleqabaGaamOAaaaaki aacIcacaWG0bGaaiykaiaadsgacaWG0baaleaacaWG0bWaa0baaWqa aiaaicdaaeaacaWGQbaaaaWcbaGaamiDamaaDaaameaacaqGMbaaba GaamOAaaaaa0Gaey4kIipakiabg2da9maalaaabaGaaGymaaqaaiaa dshadaqhaaWcbaGaaeOzaaqaaiaadQgaaaGccqGHsislcaWG0bWaa0 baaSqaaiaaicdaaeaacaWGQbaaaaaakmaapedabaGaai4waiaaigda cqGHsisldaqeWbqaaiaacIcacaaIXaGaeyOeI0IaamiuamaaDaaale aacaqGKbaabaGaamOAaaaakiaacIcacaWG0bGaaiilaiaadkhacaGG PaGaaiykaaWcbaGaamOCaiabg2da9iaaigdaaeaacaWGUbWaaSbaaW qaaiaabkhaaeqaaaqdcqGHpis1aOGaaiyxaiaadsgacaWG0baaleaa caWG0bWaa0baaWqaaiaaicdaaeaacaWGQbaaaaWcbaGaamiDamaaDa aameaacaqGMbaabaGaamOAaaaaa0Gaey4kIipakiaacYcaaaa@7CC3@           (11)

J dest j = f dest ( x j ( t f j ))=1 i n s D i j ( x j ( t f j )) / n s , MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaGaamOsamaaDa aaleaacaqGKbGaaeyzaiaabohacaqG0baabaGaamOAaaaakabbaaaa aaaa6QDYRcWdbiabg2da9iaadAgadaWgaaWcbaacbaGaa8hzaiaa=v gacaWFZbGaa8hDaaqabaGccaGGOaacbmGaa4hEamaaCaaaleqabaGa amOAaaaakiaacIcapaGaamiDamaaDaaaleaacaqGMbaabaGaamOAaa aak8qacaGGPaGaaiykaiabg2da9iaaigdacqGHsisldaaeWbqaaiaa dseadaqhaaWcbaGaamyAaaqaaiaadQgaaaGccaGGOaGaa4hEamaaCa aaleqabaGaamOAaaaakiaacIcapaGaamiDamaaDaaaleaacaqGMbaa baGaamOAaaaak8qacaGGPaGaaiykaaWcbaGaamyAaaqaaiaad6gada WgaaadbaGaam4CaaqabaaaniabggHiLdGccaGGVaGaamOBamaaBaaa leaacaWGZbaabeaakiaacYcaaaa@64BE@           (12)

where t 0 j MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaGaamiDamaaDa aaleaacaaIWaaabaGaamOAaaaaaaa@3A71@ and t f j MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaGaamiDamaaDa aaleaacaqGMbaabaGaamOAaaaaaaa@3AA0@ denote the initial and terminal time of UCAVj. The first term, J t j MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaGaamOsamaaDa aaleaacaqG0baabaGaamOAaaaaaaa@3A84@ denotes the total fight time of UCAVj Γ=L/ v min MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaGaeu4KdCKaey ypa0Jaamitaiaac+cacaWG2bWaaSbaaSqaaiGac2gacaGGPbGaaiOB aaqabaaaaa@3F8D@ . is maximum flight time along the straight line (L denotes its length) from the IP to the center of AARs. The value is predetermined, and can be used to insure the first cost term consistent with others. The second term J prd j MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaqcaaQaamOsaO Waa0baaSqaaiaabchacaqGYbGaaeizaaqaaiaadQgaaaaaaa@3D10@ , describes the average detection-probability of the nr-radar system to UCAVj, where P d j (t,r) MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaGaamiuamaaDa aaleaacaqGKbaabaGaamOAaaaakiaacIcacaWG0bGaaiilaiaadkha caGGPaaaaa@3E7E@ is the radar detection probability model between the trajectory point of UCAVj at the time t and the rth radar9. And the third term J dest j MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaqcaakeeaaaaa aaaOR2jVkapeGaamOsaOWaa0baaSqaaGqaaiaa=rgacaWFLbGaa83C aiaa=rhaaeaapaGaamOAaaaaaaa@4138@ , is the minimal form of the target damage probability, where D i j ( x j ( t f j )) MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaaeeaaaaaaaaO R2jVkapeGaamiramaaDaaaleaacaWGPbaabaWdaiaadQgaaaGcpeGa aiikaGqadiaa=HhadaahaaWcbeqaa8aacaWGQbaaaOWdbiaacIcapa GaamiDamaaDaaaleaacaqGMbaabaGaamOAaaaak8qacaGGPaGaaiyk aaaa@45ED@ denotes the target damage probability of UCAVj from the ith Monte Carlo simulation, and ns is the total simulation times.

2.4   Cooperative Trajectory Planning Problem Formulation


After establishing the above three models, the cooperative trajectory planning problem is formulated in the subsection. The problem under consideration is a cooperative scenario, consisting of Nv similar UCAVs, and the dynamics of UCAVj is given by

x ˙ j (t)= f j [ x j (t), u j (t),t], ( j=1,, N v ), MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaacbmGab8hEay aacaWaaWbaaSqabeaacaWGQbaaaOGaaiikaiaadshacaGGPaGaeyyp a0JaamOzamaaCaaaleqabaGaamOAaaaakiaacUfacaWF4bWaaWbaaS qabeaacaWGQbaaaOGaaiikaiaadshacaGGPaGaaiilaiaa=vhadaah aaWcbeqaaiaadQgaaaGccaGGOaGaamiDaiaacMcacaGGSaGaamiDai aac2facaGGSaGaaeiiaiaabIcacaqGGaGaamOAaiaab2dacaqGXaGa aeilaiablAciljaacYcacaWGobWaaSbaaSqaaiaabAhaaeqaaOGaai ykaiaacYcaaaa@571B@     (13)

where x= [x,y,h,v,γ,ψ] T 6 MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaGaaCiEaiabg2 da9iaacUfacaWG4bGaaiilaiaadMhacaGGSaGaamiAaiaacYcacaWG 2bGaaiilaiabeo7aNjaacYcacqaHipqEcaGGDbWaaWbaaSqabeaaca qGubaaaOGaeyicI4SaeSyhHe6aaWbaaSqabeaacaaI2aaaaaaa@4B20@ and u= [μ, n x , n h ] T 3 MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaGaaCyDaiabg2 da9iaacUfacqaH8oqBcaGGSaGaamOBamaaBaaaleaacaqG4baabeaa kiaacYcacaWGUbWaaSbaaSqaaiaabIgaaeqaaOGaaiyxamaaCaaale qabaGaaeivaaaakiabgIGiolabl2riHoaaCaaaleqabaGaaG4maaaa aaa@47A0@ denote the state and control vectors, which are in accordance with the following UCAV kinematic and dynamics equations10

x ˙ =vcosγcosψ,  y ˙ =vcosγsinψ,  h ˙ =vsinγ, MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbqaG=lqadIhaga Gaaiabg2da9iaadAhaciGGJbGaai4BaiaacohacqaHZoWzciGGJbGa ai4BaiaacohacqaHipqEcaGGSaGaaeiiaiqadMhagaGaaiabg2da9i aadAhaciGGJbGaai4BaiaacohacqaHZoWzciGGZbGaaiyAaiaac6ga cqaHipqEcaGGSaGaaeiiaiqadIgagaGaaiabg2da9iaadAhaciGGZb GaaiyAaiaac6gacqaHZoWzcaGGSaaaaa@5C3B@     (14)

v ˙ =g( n x sinγ),  γ ˙ = g v ( n h cosμcosγ),  ψ ˙ = g vcosγ n h sinμ, MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaGabmODayaaca Gaeyypa0Jaam4zaiaacIcacaWGUbWaaSbaaSqaaiaabIhaaeqaaOGa eyOeI0Iaci4CaiaacMgacaGGUbGaeq4SdCMaaiykaiaacYcacaqGGa Gafq4SdCMbaiaacqGH9aqpdaWcaaqaaiaadEgaaeaacaWG2baaaiaa cIcacaWGUbWaaSbaaSqaaiaabIgaaeqaaOGaci4yaiaac+gacaGGZb GaeqiVd0MaeyOeI0Iaci4yaiaac+gacaGGZbGaeq4SdCMaaiykaiaa cYcacaqGGaGafqiYdKNbaiaacqGH9aqpdaWcaaqaaiaadEgaaeaaca WG2bGaci4yaiaac+gacaGGZbGaeq4SdCgaaiaad6gadaWgaaWcbaGa aeiAaaqabaGcciGGZbGaaiyAaiaac6gacqaH8oqBcaGGSaaaaa@68D1@      (15)

where x, y, h are the east, north, and up components of the earth-fixed reference frame, and denote longitude, latitude and altitude respectively; g denotes the gravity acceleration.Fig 1. shows the spatial relation of the states.

As mentioned above, to generate optimal or suboptimal cooperative trajectories, the trajectory planning problem for the CA/GTA missions can be formulated as a cooperative trajectory optimal control problem (CTP-OCP).


Figure 1. Relation of aircraft position, velocity, flight-path and heading angels.



2.4.1 Problem 1 (CTP-OCP)

Find the trajectories, which drive the system from given initial conditions to desired final conditions over time horizons [t0, tf], while the cooperative objective function is minimized as follows

minJ= j=1 N v J j (x,u) = j=1 N v [ w t j ( t f j t 0 j )/Γ+ w r j t 0 j t f j PRD j (t)dt /( t f j t 0 j )+ w d j f dest ( x j ( t f j ) )] MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGceaqabeaaciGGTb GaaiyAaiaac6gacaWGkbGaeyypa0ZaaabCaeaacaWGkbWaaWbaaSqa beaacaWGQbaaaOGaaiikaiaahIhacaGGSaGaaCyDaiaacMcaaSqaai aadQgacqGH9aqpcaaIXaaabaGaamOtamaaBaaameaacaqG2baabeaa a0GaeyyeIuoakiabg2da9maaqahabaGaai4waiaadEhadaqhaaWcba GaaeiDaaqaaiaadQgaaaGccaGGOaGaamiDamaaDaaaleaacaqGMbaa baGaamOAaaaakiabgkHiTiaadshadaqhaaWcbaGaaGimaaqaaiaadQ gaaaGccaGGPaGaai4laiabfo5ahjabgUcaRiaadEhadaqhaaWcbaGa aeOCaaqaaiaadQgaaaGcdaWdXaqaaiaabcfacaqGsbGaaeiramaaCa aaleqabaGaamOAaaaakiaacIcacaWG0bGaaiykaiaadsgacaWG0baa leaacaWG0bWaa0baaWqaaiaaicdaaeaacaWGQbaaaaWcbaGaamiDam aaDaaameaacaqGMbaabaGaamOAaaaaa0Gaey4kIipakiaac+cacaGG OaGaamiDamaaDaaaleaacaqGMbaabaGaamOAaaaakiabgkHiTiaads hadaqhaaWcbaGaaGimaaqaaiaadQgaaaGccaGGPaGaey4kaScaleaa caWGQbGaeyypa0JaaGymaaqaaiaad6eadaWgaaadbaGaaeODaaqaba aaniabggHiLdaakeaacaaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaa ykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaG PaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaM c8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaayk W7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPa VlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8 UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7 caaMc8UaaGPaVlaadEhadaqhaaWcbaGaaeizaaqaaiaadQgaaaGcqq aaaaaaaaGUAN8Qa8qacaWGMbWaaSbaaSqaaGqaaiaa=rgacaWFLbGa a83Caiaa=rhaaeqaaOWaaeWaaeaaieWacaGF4bWdamaaCaaaleqaba GaamOAaaaak8qacaGGOaWdaiaadshadaqhaaWcbaGaaeOzaaqaaiaa dQgaaaGcpeGaaiykaaGaayjkaiaawMcaaiaac2faaaaa@E7B0@      (16)

subject to the dynamics equation (i.e., Eqn (13)) and the boundary constraints (i.e., the initial and terminal states (i.e., Eqn (6))

Φ 0 [ x j ( t 0 j ), u j ( t 0 j ), t 0 j ]=0, Φ 1 [ x j ( t f j ), u j ( t f j ), t f j ]0, MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGceaqabeaacqqHMo GrdaWgaaWcbaGaaGimaaqabaGccaGGBbacbmGaa8hEamaaCaaaleqa baGaamOAaaaakiaacIcacaWG0bWaa0baaSqaaiaaicdaaeaacaWGQb aaaOGaaiykaiaacYcacaWF1bWaaWbaaSqabeaacaWGQbaaaOGaaiik aiaadshadaqhaaWcbaGaaGimaaqaaiaadQgaaaGccaGGPaGaaiilai aadshadaqhaaWcbaGaaGimaaqaaiaadQgaaaGccaGGDbGaeyypa0Ja aGimaiaacYcaaeaacqqHMoGrdaWgaaWcbaGaaGymaaqabaGccaGGBb Gaa8hEamaaCaaaleqabaGaamOAaaaakiaacIcacaWG0bWaa0baaSqa aiaabAgaaeaacaWGQbaaaOGaaiykaiaacYcacaWF1bWaaWbaaSqabe aacaWGQbaaaOGaaiikaiaadshadaqhaaWcbaGaaeOzaaqaaiaadQga aaGccaGGPaGaaiilaiaadshadaqhaaWcbaGaaeOzaaqaaiaadQgaaa GccaGGDbGaeyizImQaaeimaiaabYcaaaaa@67CE@      (17)

and several inequality and equality constraints, the individual and cooperative constraints, including the state and control (i.e., Eqns (3)-(5) and Eqns (7)-(8)), are denoted as:



C 1 j [ x j (t), u j (t),t]0, C 2 j [ x j (t), u j (t), x k (t), u k (t),t]0 , f T j ( t f j , t f k )0, ( j,k=1,2,, N v ,jk). MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGceaqabeaacaWGdb Waa0baaSqaaiaaigdaaeaacaWGQbaaaOGaai4waGqadiaa=Hhadaah aaWcbeqaaiaadQgaaaGccaGGOaGaamiDaiaacMcacaGGSaGaa8xDam aaCaaaleqabaGaamOAaaaakiaacIcacaWG0bGaaiykaiaacYcacaWG 0bGaaiyxaiabgsMiJkaaicdacaGGSaaabaGaam4qamaaDaaaleaaca aIYaaabaGaamOAaaaakiaacUfacaWF4bWaaWbaaSqabeaacaWGQbaa aOGaaiikaiaadshacaGGPaGaaiilaiaa=vhadaahaaWcbeqaaiaadQ gaaaGccaGGOaGaamiDaiaacMcacaGGSaGaa8hEamaaCaaaleqabaGa am4AaaaakiaacIcacaWG0bGaaiykaiaacYcacaWF1bWaaWbaaSqabe aacaWGRbaaaOGaaiikaiaadshacaGGPaGaaiilaiaadshacaGGDbGa eyizImQaaGimaiaabccacaGGSaaabaGaamOzamaaDaaaleaacaqGub aabaGaamOAaaaakiaacIcacaWG0bWaa0baaSqaaiaabAgaaeaacaWG QbaaaOGaaiilaiaadshadaqhaaWcbaGaaeOzaaqaaiaadUgaaaGcca GGPaGaeyizImQaaGimaiaacYcacaqGGaGaaeikaiaabccacaWGQbGa aiilaiaadUgacqGH9aqpcaaIXaGaaiilaiaaikdacaGGSaGaeS47IW Kaaiilaiaad6eadaWgaaWcbaGaaeODaaqabaGccaGGSaGaamOAaiab gcMi5kaadUgacaGGPaGaaiOlaaaaaa@88EE@      (18)

Note that, the terminal time tf is free. Thus the problem is a free terminal time problem. In addition, it is important in the problem formulation to scale the variables. The choice of scaling will balance the equations for numerical analysis, thus improving the accuracy of the solution and the computation time11.

To solve the previous CTP-OCP as quickly as possible, an efficient numerically cooperative trajectory planning algorithm is introduced, which combines several classical techniques, including differential flatness theory, B-spline curves, and nonlinear programming. In addition, for the time-critical cooperative missions, a novel strategy is proposed to handle the temporal constraints.

3.1  Time Cooperative Strategy

The time factor of trajectories is an argument of the state and control. To deal with temporal constraints, the time along the trajectories should be considered separately. In the work, the independent intermediate variable (called virtual time here τ[0,1] MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaGaeqiXdqNaey icI4Saai4waiaaicdacaGGSaGaaGymaiaac2faaaa@3ED0@ ) is introduced and described as


τ (t t 0 )/ ( t f t 0 ) . MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaGaeqiXdqNaeS ixIa0aaSGbaeaacaGGOaGaamiDaiabgkHiTiaadshadaWgaaWcbaGa aGimaaqabaGccaGGPaaabaGaaiikaiaadshadaWgaaWcbaGaaeOzaa qabaGccqGHsislcaWG0bWaaSbaaSqaaiaaicdaaeqaaOGaaiykaaaa caGGUaaaaa@46ED@       (19)

Such that the trajectories can be generated in the virtual time domain from τ 0 =0 MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaGaeqiXdq3aaS baaSqaaiaaicdaaeqaaOGaeyypa0JaaGimaaaa@3C17@ to τ f =1 MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaGaeqiXdq3aaS baaSqaaiaabAgaaeqaaOGaeyypa0JaaGymaaaa@3C47@ .


For UCAVs are required to takeoff at the same time, it is assumed that the initial time of all UCAVs is zero (i.e. t 0 =0 MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaGaamiDamaaBa aaleaacaqGWaaabeaakiabg2da9iaaicdaaaa@3B44@ ). Thus, the terminal time t f MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipy0lf9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabeqaamaabaabaaGcbaGaamiDamaaBa aaleaacaqGMbaabeaaaaa@3786@ can be written as t f =t/τ MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaGaamiDamaaBa aaleaacaqGMbaabeaakiabg2da9maalyaabaGaamiDaaqaaiabes8a 0baaaaa@3D94@ . That is t f MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipy0lf9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabeqaamaabaabaaGcbaGaamiDamaaBa aaleaacaqGMbaabeaaaaa@3786@ , denotes the ratio between the true time variable and the newly defined virtual time variable. To coordinate the arrival time of all UCAVs, the terminal time can be defined as an argument in the dynamics to be optimized, designated as T f MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaGaamivamaaBa aaleaacaqGMbaabeaaaaa@3990@ . Then, the following relationship between the virtual time and true time domain can be obtained for an arbitrary variable χ MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaGaeq4Xdmgaaa@3959@



χ ˙ (t)= dχ(t) dt = dχ(t) T f dτ = χ (τ)/ T f . MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaGafq4XdmMbai aacaGGOaGaamiDaiaacMcacqGH9aqpdaWcaaqaaiaadsgacqaHhpWy caGGOaGaamiDaiaacMcaaeaacaWGKbGaamiDaaaacqGH9aqpdaWcaa qaaiaadsgacqaHhpWycaGGOaGaamiDaiaacMcaaeaacaWGubWaaSba aSqaaiaabAgaaeqaaOGaamizaiabes8a0baacqGH9aqpdaWcgaqaai qbeE8aJzaafaGaaiikaiabes8a0jaacMcaaeaacaWGubWaaSbaaSqa aiaabAgaaeqaaaaakiaac6caaaa@56F3@        (20)

Especially, the derivative of the speed variable can be denoted as

v(t)= x ˙ (t) 2 + y ˙ (t) 2 + h ˙ (t) 2 = ( x ) 2 + ( y ) 2 + ( h ) 2 / T f = v(τ)/ T f MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaGaamODaiaacI cacaWG0bGaaiykaiabg2da9maakaaabaGabmiEayaacaGaaiikaiaa dshacaGGPaWaaWbaaSqabeaacaaIYaaaaOGaey4kaSIabmyEayaaca GaaiikaiaadshacaGGPaWaaWbaaSqabeaacaaIYaaaaOGaey4kaSIa bmiAayaacaGaaiikaiaadshacaGGPaWaaWbaaSqabeaacaaIYaaaaa qabaGccqGH9aqpdaWcgaqaamaakaaabaGaaiikaiqadIhagaqbaiaa cMcadaahaaWcbeqaaiaaikdaaaGccqGHRaWkcaGGOaGabmyEayaafa GaaiykamaaCaaaleqabaGaaGOmaaaakiabgUcaRiaacIcaceWGObGb auaacaGGPaWaaWbaaSqabeaacaaIYaaaaaqabaaakeaacaWGubWaaS baaSqaaiaabAgaaeqaaaaakiabg2da9maalyaabaGaamODaiaacIca cqaHepaDcaGGPaaabaGaamivamaaBaaaleaacaqGMbaabeaaaaaaaa@6096@        (21)

where the superscript () represents the derivative with respect to the virtual time.
According to Eqns (13) and (19), the dynamics can be rewritten as

x j (τ)= T f x ˙ j (t)= T f f j [ x j (t), u j (t)]= f ^ j [ x j (τ), u j (τ), T f j ] T f j =0,          ( j=1,, N v ) MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGceaabbeaaceWH4b GbauaadaahaaWcbeqaaiaadQgaaaGccaGGOaGaeqiXdqNaaiykaiab g2da9iaadsfadaWgaaWcbaGaaeOzaaqabaGccqGHflY1ieWaceWF4b GbaiaadaahaaWcbeqaaiaadQgaaaGccaGGOaGaamiDaiaacMcacqGH 9aqpcaWGubWaaSbaaSqaaiaabAgaaeqaaOGaeyyXICTaamOzamaaCa aaleqabaGaamOAaaaakiaacUfacaWF4bWaaWbaaSqabeaacaWGQbaa aOGaaiikaiaadshacaGGPaGaaiilaiaa=vhadaahaaWcbeqaaiaadQ gaaaGccaGGOaGaamiDaiaacMcacaGGDbGaeyypa0JabmOzayaajaWa aWbaaSqabeaacaWGQbaaaOGaai4waiaa=HhadaahaaWcbeqaaiaadQ gaaaGccaGGOaGaeqiXdqNaaiykaiaacYcacaWF1bWaaWbaaSqabeaa caWGQbaaaOGaaiikaiabes8a0jaacMcacaGGSaGaamivamaaDaaale aacaqGMbaabaGaamOAaaaakiaac2faaeaacaWGubWaaSbaaSqaaiaa bAgaaeqaaOWaaWbaaSqabeaakiadacTHYaIOaaWaaWbaaSqabeaaca WGQbaaaOGaaeypaiaabcdacaqGSaGaaeiiaiaabccacaqGGaGaaeii aiaabccacaqGGaGaaeiiaiaabccacaqGGaGaaeiiaiaabIcacaqGGa GaamOAaiaab2dacaqGXaGaaeilaiablAciljaacYcacaWGobWaaSba aSqaaiaabAhaaeqaaOGaaiykaaaaaa@84A7@        (22)

Hence, the Problem 1 can be reformulated in virtual time domain (VTD) as Problem 2.

3.1.1 Problem 2 (CTP-OCP-VTD)

Minimize the cooperative cost function (Eqn (16)) of all UCAVs represented with respect to the new independent variable τ MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipy0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabmqaamaabaabaaGcbaGaeqiXdqhaaa@379E@ as:

J= j=1 N v J j (x,u, T f ) = j=1 N v [ w t j T f j /Γ+ w r j 0 1 PRD j (τ)dτ / T f j + w d j f dest ( x j (1) )] MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaGaamOsaiabg2 da9maaqahabaGaamOsamaaCaaaleqabaGaamOAaaaakiaacIcacaWH 4bGaaiilaiaahwhacaGGSaGaamivamaaBaaaleaacaqGMbaabeaaki aacMcaaSqaaiaadQgacqGH9aqpcaaIXaaabaGaamOtamaaBaaameaa caqG2baabeaaa0GaeyyeIuoakiabg2da9maaqahabaGaai4waiaadE hadaqhaaWcbaGaaeiDaaqaaiaadQgaaaGccaWGubWaa0baaSqaaiaa bAgaaeaacaWGQbaaaOGaai4laiabfo5ahjabgUcaRiaadEhadaqhaa WcbaGaaeOCaaqaaiaadQgaaaGcdaWdXaqaaiaabcfacaqGsbGaaeir amaaCaaaleqabaGaamOAaaaakiaacIcacqaHepaDcaGGPaGaamizai abes8a0bWcbaGaaGimaaqaaiaaigdaa0Gaey4kIipakiaac+cacaWG ubWaa0baaSqaaiaabAgaaeaacaWGQbaaaOGaey4kaSIaam4DamaaDa aaleaacaqGKbaabaGaamOAaaaakabbaaaaaaaa6QDYRcWdbiaadAga daWgaaWcbaacbaGaa8hzaiaa=vgacaWFZbGaa8hDaaqabaGcdaqada qaaGqadiaa+HhapaWaaWbaaSqabeaacaWGQbaaaOWdbiaacIcapaGa aGyma8qacaGGPaaacaGLOaGaayzkaaWdaiaac2faaSqaaiaadQgacq GH9aqpcaaIXaaabaGaamOtamaaBaaameaacaqG2baabeaaa0Gaeyye Iuoaaaa@8184@        (23)

subject to the dynamics in Eqn (22), and the boundary constraints in Eqn (17), written as:

Φ ˜ 0 [ x j (0), u j (0)]=0 Φ ˜ 1 [ x j (1), u j (1), T f j ]0 MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGceaqabeaacuqHMo GrgaacamaaBaaaleaacaaIWaaabeaakiaacUfaieWacaWF4bWaaWba aSqabeaacaWGQbaaaOGaaiikaiaaicdacaGGPaGaaiilaiaa=vhada ahaaWcbeqaaiaadQgaaaGccaGGOaGaaGimaiaacMcacaGGDbGaeyyp a0JaaGimaaqaaiqbfA6agzaaiaWaaSbaaSqaaiaaigdaaeqaaOGaai 4waiaa=HhadaahaaWcbeqaaiaadQgaaaGccaGGOaGaaGymaiaacMca caGGSaGaa8xDamaaCaaaleqabaGaamOAaaaakiaacIcacaaIXaGaai ykaiaacYcacaWGubWaa0baaSqaaiaabAgaaeaacaWGQbaaaOGaaiyx aiabgsMiJkaabcdaaaaa@5A0C@        (24)

and the inequality and equality constraints (Eqn (18)), and additional temporal constraints

C ˜ 1 j [ x j (τ), u j (τ), T f j ]0 C ˜ 2 j [ x j (τ), u j (τ), x k (τ), u k (τ), T f j , T f k ]0 f ˜ T j ( T f j , T f k )0,  (j,k=1,2,, N v ,jk) MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGceaqabeaaceWGdb GbaGaadaqhaaWcbaGaaGymaaqaaiaadQgaaaGccaGGBbacbmGaa8hE amaaCaaaleqabaGaamOAaaaakiaacIcacqaHepaDcaGGPaGaaiilai aa=vhadaahaaWcbeqaaiaadQgaaaGccaGGOaGaeqiXdqNaaiykaiaa cYcacaWGubWaa0baaSqaaiaabAgaaeaacaWGQbaaaOGaaiyxaiabgs MiJkaaicdaaeaaceWGdbGbaGaadaqhaaWcbaGaaGOmaaqaaiaadQga aaGccaGGBbGaa8hEamaaCaaaleqabaGaamOAaaaakiaacIcacqaHep aDcaGGPaGaaiilaiaa=vhadaahaaWcbeqaaiaadQgaaaGccaGGOaGa eqiXdqNaaiykaiaacYcacaWF4bWaaWbaaSqabeaacaWGRbaaaOGaai ikaiabes8a0jaacMcacaGGSaGaa8xDamaaCaaaleqabaGaam4Aaaaa kiaacIcacqaHepaDcaGGPaGaaiilaiaadsfadaqhaaWcbaGaaeOzaa qaaiaadQgaaaGccaGGSaGaamivamaaDaaaleaacaqGMbaabaGaam4A aaaakiaac2facqGHKjYOcaaIWaaabaGabmOzayaaiaWaa0baaSqaai aabsfaaeaacaWGQbaaaOGaaiikaiaadsfadaqhaaWcbaGaaeOzaaqa aiaadQgaaaGccaGGSaGaamivamaaDaaaleaacaqGMbaabaGaam4Aaa aakiaacMcacqGHKjYOcaaIWaGaaiilaiaabccacaqGGaGaaeikaiaa dQgacaGGSaGaam4Aaiabg2da9iaaigdacaGGSaGaaGOmaiaacYcacq WIVlctcaGGSaGaamOtamaaBaaaleaacaqG2baabeaakiaacYcacaWG QbGaeyiyIKRaam4AaiaacMcaaaaa@9265@        (25)

3.2  Problem Formulation in the Output Space

Due to the complexity in solving this higher dimensional space system, the differential flatness approach is introduced to transform the system dynamics to a lower dimensional space12. To confirm that the dynamics system is differentially flat, the spatial trajectory ρ(τ) MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaGaeqyWdiNaai ikaiabes8a0jaacMcaaaa@3C80@ in virtual time domain and the terminal time T f MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaciaabmqaamaabaabaaGcbaGaamivamaaBa aaleaacaqGMbaabeaaaaa@37C8@ as the flat output vector can be defined as


z=[ρ(τ), T f ]= [x(τ),y(τ),h(τ), T f ] T MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaacbaGaa8NEai abg2da9iaacUfacqaHbpGCcaGGOaGaeqiXdqNaaiykaiaacYcacaWG ubWaaSbaaSqaaiaabAgaaeqaaOGaaiyxaiabg2da9iaacUfacaWG4b Gaaiikaiabes8a0jaacMcacaGGSaGaamyEaiaacIcacqaHepaDcaGG PaGaaiilaiaadIgacaGGOaGaeqiXdqNaaiykaiaacYcacaWGubWaaS baaSqaaiaabAgaaeqaaOGaaiyxamaaCaaaleqabaGaamivaaaaaaa@5708@        (26)

The original state x and control u can be recovered from the flat outputs and their derivatives as follows13

ξ=(z, z ),    x=φ(ξ),   u=α(ξ) MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbiqaaG4ccqaH+o aEcqGH9aqpcaGGOaGaaCOEaiaacYcaceWH6bGbauaacaGGPaGaaiil aiaabccacaqGGaGaaeiiaiaabccaieWacaWF4bGaeyypa0JaeqOXdO Maaiikaiabe67a4jaacMcacaGGSaGaaeiiaiaabccacaqGGaGaa8xD aiabg2da9iabeg7aHjaacIcacqaH+oaEcaGGPaaaaa@528C@        (27)

According to the kinematics equations in Eqn (14), the remaining three states using the flatness outputs in virtual time domain can be easily described as below

v= ( x ) 2 + ( y ) 2 + ( h ) 2  , ψ=arctan( y / x ), γ=arctan( h / ( x ) 2 + ( y ) 2 ) MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaGaamODaiabg2 da9maakaaabaGaaiikaiqadIhagaqbaiaacMcadaahaaWcbeqaaiaa ikdaaaGccqGHRaWkcaGGOaGabmyEayaafaGaaiykamaaCaaaleqaba GaaGOmaaaakiabgUcaRiaacIcaceWGObGbauaacaGGPaWaaWbaaSqa beaacaaIYaaaaaqabaGccaqGGaGaaiilaiaabccacqaHipqEcqGH9a qpciGGHbGaaiOCaiaacogacaGG0bGaaiyyaiaac6gadaqadaqaamaa lyaabaGabmyEayaafaaabaGabmiEayaafaaaaaGaayjkaiaawMcaai aacYcacaqGGaGaeq4SdCMaeyypa0JaciyyaiaackhacaGGJbGaaiiD aiaacggacaGGUbWaaeWaaeaadaWcgaqaaiqadIgagaqbaaqaamaaka aabaGaaiikaiqadIhagaqbaiaacMcadaahaaWcbeqaaiaaikdaaaGc cqGHRaWkcaGGOaGabmyEayaafaGaaiykamaaCaaaleqabaGaaGOmaa aaaeqaaaaaaOGaayjkaiaawMcaaaaa@673E@        (28)

Thus, according to Eqn (15), the control variables can be determined as

n x = v / ( T f 2 g) +sinγ,   n h = ( v γ + T f 2 gcosγ ) 2 + ( v ψ cosγ ) 2 / T f g , μ=arctan( v ψ cosγ/ (v γ + T f 2 gcosγ) ), MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGceaabbeaacaWGUb WaaSbaaSqaaiaabIhaaeqaaOGaeyypa0ZaaSGbaeaaceWG2bGbauaa aeaacaGGOaGaamivamaaDaaaleaacaqGMbaabaGaaGOmaaaakiaadE gacaGGPaaaaiabgUcaRiGacohacaGGPbGaaiOBaiabeo7aNjaacYca caqGGaGaaeiiaiaad6gadaWgaaWcbaGaaeiAaaqabaGccqGH9aqpda WcgaqaamaakaaabaWaaeWaaeaacaWG2bGafq4SdCMbauaacqGHRaWk caWGubWaa0baaSqaaiaabAgaaeaacaaIYaaaaOGaam4zaiGacogaca GGVbGaai4Caiabeo7aNbGaayjkaiaawMcaamaaCaaaleqabaGaaGOm aaaakiabgUcaRmaabmaabaGaamODaiqbeI8a5zaafaGaci4yaiaac+ gacaGGZbGaeq4SdCgacaGLOaGaayzkaaWaaWbaaSqabeaacaaIYaaa aaqabaaakeaacaWGubWaaSbaaSqaaiaabAgaaeqaaOGaam4zaaaaca qGSaaabaGaeqiVd0Maeyypa0JaciyyaiaackhacaGGJbGaaiiDaiaa cggacaGGUbWaaeWaaeaadaWcgaqaaiaadAhacuaHipqEgaqbaiGaco gacaGGVbGaai4Caiabeo7aNbqaaiaacIcacaWG2bGafq4SdCMbauaa cqGHRaWkcaWGubWaa0baaSqaaiaabAgaaeaacaaIYaaaaOGaam4zai GacogacaGGVbGaai4Caiabeo7aNjaacMcaaaaacaGLOaGaayzkaaGa aiilaaaaaa@86DC@        (29)

Obviously, the dynamics constraints of this system (i.e. Eqns (14-15)) can be automatically satisfied. And the system of cooperative planning of UCAVs, including the objective function and the constraints, is mapped to a lower dimensional output space. Hence, Problem 2 can be modified as problem 3.


3.2.1 Problem 3 (CTP-OCP-VTD in Output Space)
Solving the problem

minJ= j=1 N v J j (ξ(τ))= j=1 N v [ ϕ j (φ( ξ j (0)),α( ξ j (0)),φ( ξ j (1)),α( ξ j (1)), T f j ) + 0 1 L j (φ( ξ j (τ)),α( ξ j (τ)), T f j )dτ )] MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGceaqabeaaciGGTb GaaiyAaiaac6gacaWGkbGaeyypa0ZaaabCaeaacaWGkbWaaWbaaSqa beaacaWGQbaaaaqaaiaadQgacqGH9aqpcaaIXaaabaGaamOtamaaBa aameaacaqG2baabeaaa0GaeyyeIuoakiaacIcacqaH+oaEcaGGOaGa eqiXdqNaaiykaiaacMcacqGH9aqpdaaeWbqaaiaacUfacqaHvpGzda ahaaWcbeqaaiaadQgaaaGccaGGOaGaeqOXdOMaaiikaiabe67a4naa CaaaleqabaGaamOAaaaakiaacIcacaaIWaGaaiykaiaacMcacaGGSa GaeqySdeMaaiikaiabe67a4naaCaaaleqabaGaamOAaaaakiaacIca caaIWaGaaiykaiaacMcacaGGSaGaeqOXdOMaaiikaiabe67a4naaCa aaleqabaGaamOAaaaakiaacIcacaaIXaGaaiykaiaacMcacaGGSaGa eqySdeMaaiikaiabe67a4naaCaaaleqabaGaamOAaaaakiaacIcaca aIXaGaaiykaiaacMcacaGGSaGaamivamaaDaaaleaacaqGMbaabaGa amOAaaaakiaacMcaaSqaaiaadQgacqGH9aqpcaaIXaaabaGaamOtam aaBaaameaacaqG2baabeaaa0GaeyyeIuoakiabgUcaRaqaaiaaykW7 caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVl aaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8Ua aGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7ca aMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaa ykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaG PaVlaaykW7caaMc8UaaGPaVlaaykW7caaMc8UaaGPaVlaaykW7caaM c8UaaGPaVlaaykW7caaMc8UaaGPaVpaapedabaGaamitamaaCaaale qabaGaamOAaaaakiaacIcacqaHgpGAcaGGOaGaeqOVdG3aaWbaaSqa beaacaWGQbaaaOGaaiikaiabes8a0jaacMcacaGGPaGaaiilaiabeg 7aHjaacIcacqaH+oaEdaahaaWcbeqaaiaadQgaaaGccaGGOaGaeqiX dqNaaiykaiaacMcacaGGSaGaamivamaaDaaaleaacaqGMbaabaGaam OAaaaakiaacMcacaWGKbGaeqiXdqhaleaacaaIWaaabaGaaGymaaqd cqGHRiI8aOGaaiykaiaac2faaaaa@F402@        (30)

subject to

L 0 j c 0 j (φ( ξ j (0)),α( ξ j (0)), T f j ) U 0 j L f j c f j (φ( ξ j (1)),α( ξ j (1)), T f j ) U f j L τ j c τ j (φ( ξ j (τ)),α( ξ j (τ)), T f j ) U τ j ,  0<τ<1 ϖ j [φ( ξ j (τ)),α( ξ j (τ)),φ( ξ k (τ)),α( ξ k (τ)), T f j , T f k ]0, f ˜ T j ( T f j , T f k )0,  (j,k=1,2,, N v ,jk). MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGceaqabeaacaWGmb Waa0baaSqaaiaaicdaaeaacaWGQbaaaOGaeyizImQaam4yamaaDaaa leaacaaIWaaabaGaamOAaaaakiaacIcacqaHgpGAcaGGOaGaeqOVdG 3aaWbaaSqabeaacaWGQbaaaOGaaiikaiaaicdacaGGPaGaaiykaiaa cYcacqaHXoqycaGGOaGaeqOVdG3aaWbaaSqabeaacaWGQbaaaOGaai ikaiaaicdacaGGPaGaaiykaiaacYcacaWGubWaa0baaSqaaiaabAga aeaacaWGQbaaaOGaaiykaiabgsMiJkaadwfadaqhaaWcbaGaaGimaa qaaiaadQgaaaaakeaacaWGmbWaa0baaSqaaiaabAgaaeaacaWGQbaa aOGaeyizImQaam4yamaaDaaaleaacaqGMbaabaGaamOAaaaakiaacI cacqaHgpGAcaGGOaGaeqOVdG3aaWbaaSqabeaacaWGQbaaaOGaaiik aiaaigdacaGGPaGaaiykaiaacYcacqaHXoqycaGGOaGaeqOVdG3aaW baaSqabeaacaWGQbaaaOGaaiikaiaaigdacaGGPaGaaiykaiaacYca caWGubWaa0baaSqaaiaabAgaaeaacaWGQbaaaOGaaiykaiabgsMiJk aadwfadaqhaaWcbaGaaeOzaaqaaiaadQgaaaaakeaacaWGmbWaa0ba aSqaaiabes8a0bqaaiaadQgaaaGccqGHKjYOcaWGJbWaa0baaSqaai abes8a0bqaaiaadQgaaaGccaGGOaGaeqOXdOMaaiikaiabe67a4naa CaaaleqabaGaamOAaaaakiaacIcacqaHepaDcaGGPaGaaiykaiaacY cacqaHXoqycaGGOaGaeqOVdG3aaWbaaSqabeaacaWGQbaaaOGaaiik aiabes8a0jaacMcacaGGPaGaaiilaiaadsfadaqhaaWcbaGaaeOzaa qaaiaadQgaaaGccaGGPaGaeyizImQaamyvamaaDaaaleaacqaHepaD aeaacaWGQbaaaOGaaiilaiaabccacaqGGaGaaGimaiabgYda8iabes 8a0jabgYda8iaaigdaaeaacqaHwpGDdaahaaWcbeqaaiaadQgaaaGc caGGBbGaeqOXdOMaaiikaiabe67a4naaCaaaleqabaGaamOAaaaaki aacIcacqaHepaDcaGGPaGaaiykaiaacYcacqaHXoqycaGGOaGaeqOV dG3aaWbaaSqabeaacaWGQbaaaOGaaiikaiabes8a0jaacMcacaGGPa GaaiilaiabeA8aQjaacIcacqaH+oaEdaahaaWcbeqaaiaadUgaaaGc caGGOaGaeqiXdqNaaiykaiaacMcacaGGSaGaeqySdeMaaiikaiabe6 7a4naaCaaaleqabaGaam4AaaaakiaacIcacqaHepaDcaGGPaGaaiyk aiaacYcacaWGubWaa0baaSqaaiaabAgaaeaacaWGQbaaaOGaaiilai aadsfadaqhaaWcbaGaaeOzaaqaaiaadUgaaaGccaGGDbGaeyizImQa aGimaiaacYcaaeaaceWGMbGbaGaadaqhaaWcbaGaaeivaaqaaiaadQ gaaaGccaGGOaGaamivamaaDaaaleaacaqGMbaabaGaamOAaaaakiaa cYcacaWGubWaa0baaSqaaiaabAgaaeaacaWGRbaaaOGaaiykaiabgs MiJkaaicdacaGGSaGaaeiiaiaabccacaqGOaGaamOAaiaacYcacaWG RbGaeyypa0JaaGymaiaacYcacaaIYaGaaiilaiabl+UimjaacYcaca WGobWaaSbaaSqaaiaabAhaaeqaaOGaaiilaiaadQgacqGHGjsUcaWG RbGaaiykaiaac6caaaaa@FE71@        (31)

3.3  Parameterization of the Spatial Trajectory

To generate feasible trajectories in a finite parameter space, flatness outputs are parameterized in terms of B-spline curves, which have been used to represent the trajectory of UCAV to minimize computation loads and successfully applied to path/trajectory planning14. Without loss of generality, to decrease the computational time, the 3-degree B-spline curves are chosen to represent the trajectory, whose order equals to 4. Therefore, choosing every four neighboring points from the point set as control points, a 3-degree B-spline curve can be defined, which is a trajectory segment. In addition, it is clearly seen from Eqn (19) that the non-decreasing variable τ MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipy0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabeqaamaabaabaaGcbaGaeqiXdqhaaa@379C@ can be the knot sequence of B-spline. Consequently, the ith trajectory segment can be described as follows


B i (τ)= b 1,4 i (τ) C 1 i + b 2,4 i (τ) C 2 i + b 3,4 i (τ) C 3 i + b 4,4 i (τ) C 4 i MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbiqaaWyacaWGcb WaaSbaaSqaaiaadMgaaeqaaOGaaiikaiabes8a0jaacMcacqGH9aqp caWGIbWaa0baaSqaaiaaigdacaGGSaGaaGinaaqaaiaadMgaaaGcca GGOaGaeqiXdqNaaiykaiaadoeadaqhaaWcbaGaaGymaaqaaiaadMga aaGccqGHRaWkcaWGIbWaa0baaSqaaiaaikdacaGGSaGaaGinaaqaai aadMgaaaGccaGGOaGaeqiXdqNaaiykaiaadoeadaqhaaWcbaGaaGOm aaqaaiaadMgaaaGccqGHRaWkcaWGIbWaa0baaSqaaiaaiodacaGGSa GaaGinaaqaaiaadMgaaaGccaGGOaGaeqiXdqNaaiykaiaadoeadaqh aaWcbaGaaG4maaqaaiaadMgaaaGccqGHRaWkcaWGIbWaa0baaSqaai aaisdacaGGSaGaaGinaaqaaiaadMgaaaGccaGGOaGaeqiXdqNaaiyk aiaadoeadaqhaaWcbaGaaGinaaqaaiaadMgaaaaaaa@68D5@        (32)

Obviously, as b m,4 i (τ) (m=1,2,3,4) MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaGaamOyamaaDa aaleaacaWGTbGaaiilaiaaisdaaeaacaWGPbaaaOGaaiikaabaaaaa aaaapeGaeqiXdq3daiaacMcacaqGGaGaaiikaiaad2gacqGH9aqpca aIXaGaaiilaiaaikdacaGGSaGaaG4maiaacYcacaaI0aGaaiykaaaa @4851@ is known, once the control point serial C m i MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaGaam4qamaaDa aaleaacaWGTbaabaGaamyAaaaaaaa@3A77@ is determined, the ith trajectory segment can be generated by Eqn (32). Thus, the trajectory of UCAV is mapped into the control point serial of B-spline curves.

As is previously stated, the spatial trajectory segment B i (τ) MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaGaamOqamaaBa aaleaacaWGPbaabeaakiaacIcaqaaaaaaaaaWdbiabes8a09aacaGG Paaaaa@3CDA@ can be replaced by flatness output variables, [ x i (τ), y i (τ), h i (τ)] T MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaGaai4waiaadI hadaWgaaWcbaGaamyAaaqabaGccaGGOaaeaaaaaaaaa8qacqaHepaD paGaaiykaiaacYcacaWG5bWaaSbaaSqaaiaadMgaaeqaaOGaaiika8 qacqaHepaDpaGaaiykaiaacYcacaWGObWaaSbaaSqaaiaadMgaaeqa aOGaaiika8qacqaHepaDpaGaaiykaiaac2fadaahaaWcbeqaaiaabs faaaaaaa@4BE1@ and simplified as

x i (τ)= a i 0 + a i 1 τ+ a i 2 τ 2 + a i 3 τ 3 ,   y i (τ)= b i 0 + b i 1 τ+ b i 2 τ 2 + b i 3 τ 3 , h i (τ)= c i 0 + c i 1 τ+ c i 2 τ 2 + c i 3 τ 3 ,  τ[0,1] MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGceaabbeaacaWG4b WaaSbaaSqaaiaadMgaaeqaaOGaaeikaabaaaaaaaaapeGaeqiXdq3d aiaabMcacqGH9aqpcaWGHbWaa0baaSqaaiaadMgaaeaacaaIWaaaaO Gaey4kaSIaamyyamaaDaaaleaacaWGPbaabaGaaGymaaaak8qacqaH epaDpaGaey4kaSIaamyyamaaDaaaleaacaWGPbaabaGaaGOmaaaak8 qacqaHepaDpaWaaWbaaSqabeaacaaIYaaaaOGaey4kaSIaamyyamaa DaaaleaacaWGPbaabaGaaG4maaaak8qacqaHepaDpaWaaWbaaSqabe aacaaIZaaaaOGaaeilaiaabccacaqGGaGaamyEamaaBaaaleaacaWG PbaabeaakiaabIcapeGaeqiXdq3daiaabMcacqGH9aqpcaWGIbWaa0 baaSqaaiaadMgaaeaacaaIWaaaaOGaey4kaSIaamOyamaaDaaaleaa caWGPbaabaGaaGymaaaak8qacqaHepaDpaGaey4kaSIaamOyamaaDa aaleaacaWGPbaabaGaaGOmaaaak8qacqaHepaDpaWaaWbaaSqabeaa caaIYaaaaOGaey4kaSIaamOyamaaDaaaleaacaWGPbaabaGaaG4maa aak8qacqaHepaDpaWaaWbaaSqabeaacaaIZaaaaOGaaeilaaqaaiaa dIgadaWgaaWcbaGaamyAaaqabaGccaqGOaWdbiabes8a09aacaqGPa Gaeyypa0Jaam4yamaaDaaaleaacaWGPbaabaGaaGimaaaakiabgUca RiaadogadaqhaaWcbaGaamyAaaqaaiaaigdaaaGcpeGaeqiXdq3dai abgUcaRiaadogadaqhaaWcbaGaamyAaaqaaiaaikdaaaGcpeGaeqiX dq3damaaCaaaleqabaGaaGOmaaaakiabgUcaRiaadogadaqhaaWcba GaamyAaaqaaiaaiodaaaGcpeGaeqiXdq3damaaCaaaleqabaGaaG4m aaaakiaabYcacaqGGaGaaeiia8qacqaHepaDcqGHiiIZcaqGBbGaae imaiaabYcacaqGXaGaaeyxaaaaaa@969E@      (33)

where a i m1 MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaGaamyyamaaDa aaleaacaWGPbaabaGaamyBaiabgkHiTiaaigdaaaaaaa@3C3D@ , b i m1 MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaGaamOyamaaDa aaleaacaWGPbaabaGaamyBaiabgkHiTiaaigdaaaaaaa@3C3E@ , c i m1 MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaGaam4yamaaDa aaleaacaWGPbaabaGaamyBaiabgkHiTiaaigdaaaaaaa@3C3F@ are the evaluation coefficients determined by the coefficients of the B-spline basis functions.

Once all control points are obtained, the spatial trajectory, i.e., the flatness outputs will be generated. Such that, the original control variables and remaining state variables, as well as the objective function and the constraints, can be parameterized (i.e., discretized) by the control point serial. Thus, the trajectory will be mapped to control point serial, and the complex trajectory planning is transformed into a parameter optimization problem.

3.4  Transformation into a Nonlinear Programming Problem

After the flatness outputs have been parameterized in terms of B-spline curves, the coefficients of the B-spline basis functions can be found by nonlinear programming.

First of all, the virtual time horizon of UCAVj denoted by τ[0,1] MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaGaeqiXdqNaey icI4Saai4waiaaicdacaGGSaGaaGymaiaac2faaaa@3ED0@ need be uniformity discretized as the ( N c j +1) MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaGaaeikaiaad6 eadaqhaaWcbaGaae4yaaqaaiaadQgaaaGccaqGRaGaaeymaiaabMca aaa@3D3A@ collocation points below

d c j =01/ N c j ... τ n j  , ... , ( N c j 1)/ N c j 1 MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaGaaeizamaaDa aaleaacaqGJbaabaGaamOAaaaakiabg2da9iaaicdacaqGSaGaaeii aiaaigdacaGGVaGaamOtamaaDaaaleaacaqGJbaabaGaamOAaaaaki aabYcacaqGGaGaaiOlaiaac6cacaGGUaGaaeilaiaabccacqaHepaD daqhaaWcbaGaamOBaaqaaiaadQgaaaGccaqGGaGaaeilaiaabccaca GGUaGaaiOlaiaac6cacaqGGaGaaeilaiaabccacaqGOaGaamOtamaa DaaaleaacaqGJbaabaGaamOAaaaakiabgkHiTiaaigdacaqGPaGaai 4laiaad6eadaqhaaWcbaGaae4yaaqaaiaadQgaaaGccaqGSaGaaeii aiaaigdaaaa@5C09@       (34)

And the coordinates of collocation points can be written as [x( τ n j ), y( τ n j ), h( τ n j )] MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaGaae4waiaadI hacaqGOaGaeqiXdq3aa0baaSqaaiaad6gaaeaacaWGQbaaaOGaaeyk aiaabYcacaqGGaGaamyEaiaabIcacqaHepaDdaqhaaWcbaGaamOBaa qaaiaadQgaaaGccaqGPaGaaeilaiaabccacaWGObGaaeikaiabes8a 0naaDaaaleaacaWGUbaabaGaamOAaaaakiaabMcacaqGDbaaaa@4E8B@ , which can be calculated by the Eqn (33). Thus, the continuous-spline curves can be discretized in ( N c j +1) MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaGaaeikaiaad6 eadaqhaaWcbaGaae4yaaqaaiaadQgaaaGccaqGRaGaaeymaiaabMca aaa@3D3A@ points, and the constraints on the B-spline curves will be enforced at collocation points.
For the 3-degree B-spline curve under consideration, B-spline coefficients of UCAVj for all position outputs can be denoted as

ζ j :=[ ( C 1 1j , C 2 1j ,, C p 1j 1j , C 1 2j , C 2 2j , , C p 2j 2j ,, C 1 qj , C 2 qj ,, C p qj qj ) ] M j MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaGaeqOTdO3aaS baaSqaaiaadQgaaeqaaOGaaiOoaiabg2da9iaacIcacqaH2oGEdaWg aaWcbaGaaGymaiaadQgaaeqaaOGaaiilaiabeA7a6naaBaaaleaaca aIYaGaamOAaaqabaGccaGGSaGaeSOjGSKaaiilaiabeA7a6naaBaaa leaacaWGXbGaamOAaaqabaGccaGGPaGaeyypa0ZaamWaaqaabeqaai aacIcacaWGdbWaa0baaSqaaiaaigdaaeaacaaIXaGaamOAaaaakiaa cYcacaWGdbWaa0baaSqaaiaaikdaaeaacaaIXaGaamOAaaaakiaacY cacqWIMaYscaGGSaGaam4qamaaDaaaleaacaWGWbWaaSbaaWqaaiaa igdacaWGQbaabeaaaSqaaiaaigdacaWGQbaaaOGaaiilaiaadoeada qhaaWcbaGaaGymaaqaaiaaikdacaWGQbaaaOGaaiilaiaadoeadaqh aaWcbaGaaGOmaaqaaiaaikdacaWGQbaaaOGaaiilaaqaaiaaykW7cq WIMaYscaGGSaGaam4qamaaDaaaleaacaWGWbWaaSbaaWqaaiaaikda caWGQbaabeaaaSqaaiaaikdacaWGQbaaaOGaaiilaiablAciljaacY cacaWGdbWaa0baaSqaaiaaigdaaeaacaWGXbGaamOAaaaakiaacYca caWGdbWaa0baaSqaaiaaikdaaeaacaWGXbGaamOAaaaakiaacYcacq WIMaYscaGGSaGaam4qamaaDaaaleaacaWGWbWaaSbaaWqaaiaadgha caWGQbaabeaaaSqaaiaadghacaWGQbaaaOGaaiykaaaacaGLBbGaay zxaaGaeyicI4SaeSyhHe6aaWbaaSqabeaacaWGnbWaaSbaaWqaaiaa dQgaaeqaaaaaaaa@896B@      (35)

where M j = i=1 q p ij MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaGaamytamaaBa aaleaacaWGQbaabeaakiabg2da9maaqahabaGaamiCamaaBaaaleaa caWGPbGaamOAaaqabaaabaGaamyAaiabg2da9iaaigdaaeaacaWGXb aaniabggHiLdaaaa@4379@ Pij is the number of the control points for the ith position output of UCAVj , and these coefficients are used as the decision variables in the NLP problem.
Thus, the multi-UCAV trajectory planning problem in optimal control framework is transformed into a NLP problem, i.e. CTP-NLP, given by Problem 4.

3.4.1 Problem 4 (CTP-NLP)
Solving the problem

min ζ K ,T N v  F(ζ,  T f ),     K= j=1 N v M j s.t.     g r (ζ,  T f )0,  r=1,2,, r n          h s (ζ,  T f )=0,  s=1,2,, s n MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGceaqabeaadaWfqa qaaiGac2gacaGGPbGaaiOBaaWcbaGaeqOTdONaeyicI4SaeSyhHe6a aWbaaWqabeaacaWGlbaaaSGaaiilaiaadsfacqGHiiIZcqWIDesOda ahaaadbeqaaiaad6eadaWgaaqaaiaabAhaaeqaaaaaaSqabaGccaqG GaGaamOraiaacIcacqaH2oGEcaGGSaGaaeiiaiaadsfadaWgaaWcba GaaeOzaaqabaGccaGGPaGaaiilaiaabccacaqGGaGaaeiiaiaabcca caqGGaGaam4saiabg2da9maaqahabaGaamytamaaBaaaleaacaWGQb aabeaaaeaacaWGQbGaeyypa0JaaGymaaqaaiaad6eadaWgaaadbaGa aeODaaqabaaaniabggHiLdaakeaacaqGZbGaaeOlaiaabshacaqGUa GaaeiiaiaabccacaqGGaGaaeiiaiaadEgadaWgaaWcbaGaamOCaaqa baGccaGGOaGaeqOTdONaaiilaiaabccacaWGubWaaSbaaSqaaiaabA gaaeqaaOGaaiykaiabgsMiJkaaicdacaGGSaGaaeiiaiaabccacaWG YbGaeyypa0JaaGymaiaacYcacaaIYaGaaiilaiablAciljaacYcaca WGYbWaaSbaaSqaaiaad6gaaeqaaaGcbaGaaeiiaiaabccacaqGGaGa aeiiaiaabccacaqGGaGaaeiiaiaabccacaWGObWaaSbaaSqaaiaado haaeqaaOGaaiikaiabeA7a6jaacYcacaqGGaGaamivamaaBaaaleaa caqGMbaabeaakiaacMcacqGH9aqpcaaIWaGaaiilaiaabccacaqGGa Gaam4Caiabg2da9iaaigdacaGGSaGaaGOmaiaacYcacqWIMaYscaGG SaGaam4CamaaBaaaleaacaWGUbaabeaaaaaa@92FF@       (36)

where ζ MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaGaeqOTdOhaaa@395F@ is a decision vector of cooperative planning, described as ζ:=( ζ 1 , ζ 2 ,, ζ N v ) MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaGaeqOTdONaai Ooaiabg2da9iaacIcacqaH2oGEdaWgaaWcbaGaaGymaaqabaGccaGG SaGaeqOTdO3aaSbaaSqaaiaaikdaaeqaaOGaaiilaiablAciljaacY cacqaH2oGEdaWgaaWcbaGaamOtamaaBaaameaacaqG2baabeaaaSqa baGccaGGPaaaaa@4902@ , r n MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaGaamOCamaaBa aaleaacaWGUbaabeaaaaa@39B8@ and s n MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaGaam4CamaaBa aaleaacaWGUbaabeaaaaa@39B9@ are the total number of the inequality and equality constraints of all UCAVs, respectively.

Then the resulting CTP-NLP problem can be solved through well developed algorithms, such as the sequential quadratic programming. In the paper, the SNOPT software toolbox is used due to its advantages in solution effectiveness for the large-scale NLP problems15.

The basic ideas presented in this paper are illustrated in the following two scenarios. The common parameters of models in our simulations are listed in Table 1. .


Table 1. State and control constraints of UCAVs



The experimental test environment is a rectangle area of 30 × 40 km2, as shown in Fig 2. and Fig. 4 All the results presented below are generated using TOMLAB/SNOPT software toolbox on a 2.4-GHz Core 2 CPU, 2G RAM computer running with MATLAB R2009b. The weighting factors, w t j MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaGaam4DamaaDa aaleaacaqG0baabaGaamOAaaaaaaa@3AB1@ , w r j MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaGaam4DamaaDa aaleaacaqGYbaabaGaamOAaaaaaaa@3AAF@ , and w d j MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaGaam4DamaaDa aaleaacaqGKbaabiqaamtbcaWGQbaaaaaa@3ADC@ are set as 0.4, 0.3, and 0.3, by using trail and error. And the minimum safety radius of each UCAV dsafe is set as 500 m. The first initial guesses of the optimization parameters can be generated by using the B-spline interpolation of the lines from the IPs to the AARs of the targets. Then, the multiresolution-based iterative strategy8 is used to generate some ‘good’ initial guesses to the solver, which can reduce the number of iterations required to solve the NLP problem. The specified maximum number of collocation points is set as. N c,max =100 MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaGaamOtamaaBa aaleaacaqGJbGaaeilaiaab2gacaqGHbGaaeiEaaqabaGccqGH9aqp caaIXaGaaGimaiaaicdaaaa@4045@ Moreover, it is assumed that the target assignment is already completed before.


Figure 2. Two collision-free UCAV trajectories with arriving simultaneously.


4.1 Scenario 1: Cooperative Trajectories of Two UCAVs Arriving Simultaneously

In this scenario, two UCAVs cooperatively attack two stationary ground targets while avoiding a series of static obstacles/threats detected and collision en route, and satisfying aircraft dynamics constraint, especially simultaneous arrival constraint. UCAV1 and UCAV2 start at each IP, i.e., IP1 (10 km, 2 km, 2 km) and IP2 (17 km, 2 km, 2 km). Then they fly into the AARs of two targets: Target 1 (4.2 km, 34 km, 0 km) and Target 2 (14 km, 40 km, 0 km), respectively. The other three initial states and control inputs are given by x ˜ 0 1 = x ˜ 0 2 =[220m/s , 0 , 90 ] MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaGabCiEayaaia Waa0baaSqaaiaaicdaaeaacaaIXaaaaOGaeyypa0JabCiEayaaiaWa a0baaSqaaiaaicdaaeaacaaIYaaaaOGaeyypa0Jaai4waiaaikdaca aIYaGaaGimaiaad2gacaGGVaGaam4CaiaabYcacaqGGaGaaeimamaa CaaaleqabaGaeSigI8gaaOGaaeilaiaabccacaqG5aGaaeimamaaCa aaleqabaGaeSigI8gaaOGaaiyxaaaa@4D5E@ and. u 0 1 = u 0 2 =[ 0 0g, 1g] MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaGaaCyDamaaDa aaleaacaaIWaaabaGaaGymaaaakiabg2da9iaahwhadaqhaaWcbaGa aGimaaqaaiaaikdaaaGccqGH9aqpcaGGBbGaaeimamaaCaaaleqaba GaeSigI8gaaOGaaeilaiaabccacaaIWaGaae4zaiaacYcacaqGGaGa aeymaiaabEgacaGGDbaaaa@48CF@ Thus Γ 1 =465 s MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaGaeu4KdC0aaW baaSqabeaacaaIXaaaaOGaeyypa0JaaGinaiaaiAdacaaI1aGaaeii aiaabohaaaa@3ED8@ and Γ 2 =530 s MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaGaeu4KdC0aaW baaSqabeaacaaIYaaaaOGaeyypa0JaaGynaiaaiodacaaIWaGaaeii aiaabohaaaa@3ED2@ . The collision-free trajectories of the UCAVs and their arrival time, i.e. the total flight time, are shown in Fig 2. In addition, the approximate weapon trajectories are drawn to simulate the attack process. The time histories of the UCAVs’ states ( (V, γ, ψ) MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaGaaiikaiaadA facaGGSaGaaeiiaiabeo7aNjaacYcacaqGGaGaeqiYdKNaaiykaaaa @3FF2@ ) and control inputs (μ,  n x ,  n h ) MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaGaaiikaiabeY 7aTjaacYcacaqGGaGaamOBamaaBaaaleaacaqG4baabeaakiaacYca caqGGaGaamOBamaaBaaaleaacaqGObaabeaakiaacMcaaaa@4190@ are shown in Fig. 3 The average detection-probabilities of UCAVs are 0.12 and 0.17, and the target damage probabilities of UCAVs are 0.851 and 0.762, respectively. It is clear that the resulting trajectories are smooth, and the constraints on these variables, especially the cooperative constraints are all satisfied (Table 1.), which means the resulting trajectories are feasible and safe.


Figure 3. State and control time histories of the two UCAVs.


4.2 Scenario 2: Cooperative Trajectories of Multi-UCAV Arriving in Sequence

In this scenario, three UCAVs attack two stationary ground targets cooperatively. The only additional requirement is that the UCAVs arrive at their AARs in sequence rather than simultaneously, and the intervals between UCAVs are equal, denoted as Δ 12 = Δ 23 =30 s MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipy0Je9sqqrpepC0xbbL8F4rqqrFfpeea0xe9Lq=Jc9 vqaqpepm0xbba9pwe9Q8fs0=yqaqpepae9pg0FirpepeKkFr0xfr=x fr=xb9adbaqaaeGaciGaaiaabeqaamaabaabaaGcbaGaeyiLdq0aaS baaSqaaiaaigdacaaIYaaabeaakiabg2da9iabgs5aenaaBaaaleaa caaIYaGaaG4maaqabaGccqGH9aqpcaaIZaGaaGimaiaabccacaqGZb aaaa@411D@ . The coordinates of IP1, IP2 and two targets are the same as Scenario 1, and the third IP is IP3 (4 km, 2 km, 2 km). The result of the previous finished target assignment is that the UCAV1 and UCAV2 attack Target 1 and UCAV3 attacks the other one. The initial remaining three states and control inputs of the UCAV3 are preset as x ˜ 0 3 =[220m/s , 0 , 90 ] MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaGabCiEayaaia Waa0baaSqaaiaaicdaaeaacaaIZaaaaOGaeyypa0Jaai4waiaaikda caaIYaGaaGimaiaad2gacaGGVaGaam4CaiaabYcacaqGGaGaaeimam aaCaaaleqabaGaeSigI8gaaOGaaeilaiaabccacaqG5aGaaeimamaa CaaaleqabaGaeSigI8gaaOGaaiyxaaaa@499E@ and u 0 3 =( 0 , 0g, 1g) MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaGaaCyDamaaDa aaleaacaaIWaaabaGaaG4maaaakiabg2da9iaacIcacaaIWaWaaWba aSqabeaacqWIyiYBaaGccaGGSaGaaeiiaiaaicdacaqGNbGaaiilai aabccacaaIXaGaae4zaiaacMcaaaa@44C7@ . Thus Γ 1 =465 s MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaGaeu4KdC0aaW baaSqabeaacaaIXaaaaOGaeyypa0JaaGinaiaaiAdacaaI1aGaaeii aiaabohaaaa@3ED8@ , Γ 2 =438 s MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaGaeu4KdC0aaW baaSqabeaacaaIYaaaaOGaeyypa0JaaGinaiaaiodacaaI4aGaaeii aiaabohaaaa@3ED9@ and Γ 3 =576 s MathType@MTEF@5@5@+= feaagKart1ev2aaatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaGaeu4KdC0aaW baaSqabeaacaaIZaaaaOGaeyypa0JaaGynaiaaiEdacaaI2aGaaeii aiaabohaaaa@3EDD@ Fig. 4 shows the overall collision-free attack trajectories of multi-UCAV and their arrival time. It can be clearly demonstrated that the UCAVs can avoid all obstacles or threats and successfully fly into the AARs in sequence to perform weapon delivery. Fig. 5 shows the distance between each pair UCAVs. From it, one can find that the minimum distance is more than the minimum safety radius of dsafe = 500 m. Fig. 6 shows the time histories of the UCAVs’ states ( (V, γ, ψ) MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaGaaiikaiaadA facaGGSaGaaeiiaiabeo7aNjaacYcacaqGGaGaeqiYdKNaaiykaaaa @3FF2@ ) and control inputs (μ,  n x ,  n h ) MathType@MTEF@5@5@+= feaagKart1ev2aqatCvAUfeBSjuyZL2yd9gzLbvyNv2CaerbuLwBLn hiov2DGi1BTfMBaeXatLxBI9gBaerbd9wDYLwzYbItLDharqqtubsr 4rNCHbGeaGqipu0JLipgYlb91rFfpec8Eeeu0xXdbba9frFj0=OqFf ea0dXdd9vqai=hGuQ8kuc9pgc9q8qqaq=dir=f0=yqaiVgFr0xfr=x fr=xb9adbaqaaeGacaGaaiaabeqaamaabaabaaGcbaGaaiikaiabeY 7aTjaacYcacaqGGaGaamOBamaaBaaaleaacaqG4baabeaakiaacYca caqGGaGaamOBamaaBaaaleaacaqGObaabeaakiaacMcaaaa@4190@ . The average detection-probabilities of UCAVs are 0.133, 0.092 and 0.12, and the target damage probabilities of UCAVs are 0.813, 0.901 and 0.724, respectively. Obviously, the resulting trajectories are feasible and safe for all the constraints listed in Table 1. are satisfied.


Figure 4. Three collision-free UCAV trajectories with arriving in sequence.



Figure 5. Distance between each pair UCAVs.



Figure 6. State and control time histories of the three UCAVs.


This paper is devoted to exploring the cooperative collision-free trajectory planning for multi-UCAV performing the CA/GTA missions. A novel cooperative planning approach is proposed in optimal control framework, based on differential flatness, B-spline curves and nonlinear programming. To integrate weapon delivery constraints in the problem formulation, an approximate AAR model is established. Moreover, the notion of the virtual time is introduced to handle the temporal constraint. Instead of solving an OCP over a high-dimensional continuous space, the NLP problem of very low dimension has been solved successfully, over a much smaller space. The proposed approach can efficiently solve the point-to-region trajectory planning problem with integrating the spatial and temporal constraints on the trajectory level, whose validity is illustrated with some simulation results finally. Further efforts may be made to analyse some uncertain factors in the true battlefield environment and carry out the research on the real-time cooperative trajectory planning.

1.     Murray, R.M. Recent research in cooperative control of multi-vehicle systems. J. Dyn. Syst.-T. ASME, 2007, 129(5), 571-583.

2.     Hurni, M. A.; Sekhavat, P.; Karpenko, M. & Ross, I. M. A pseudospectral optimal motion planner for autonomous unmanned vehicles. In American Control Conference. Baltimore, MD, USA. 2010.[Full text via CrossRef]

3.     McLain, T.W. & Beard, R.W. Coordination variables, coordination functions, and cooperative-timing missions. J. Guid. Control Dynam., 2005, 28(1), 150-161.[Full text via CrossRef]

4.     Kaminer, I.; Yakimenko, O.; Dobrokhodov, V.; A. Pascoal; Hovakimyan, N.; Cao, C.; Young, A. & Patel, V. Coordinated path pollowing for time-critical missions of multiple UAVs via L1 adaptive output feedback controllers. In AIAA Guidance, Navigation and Control Conference and Exhibit. Hilton Head, South Carolina, USA. 2007, pp.1-23.[Full text via CrossRef]

5.     Bollino, K.P. & Lewis, L.R. Collision-free multi-UAV optimal path planning and cooperative control for tactical applications. In AIAA Guidance, Navigation and Control Conference and Exhibit. Honolulu, Hawaii, USA. 2008, pp.1-18.[Full text via CrossRef]

6.     Lian, F.-L. Cooperative path planning of dynamical multi-agent systems using differential flatness approach. Int. J. Control Autom., 2008, 6(3), 401-412.

7.     Kuwata, Y. & How, J.P. Cooperative distributed robust trajectory optimization using receding horizon MILP. IEEE T. Contr. Syst. T., 2011, 19(2), 423-431.[Full text via CrossRef]

8.     Zhang, Y.; Chen, J. & Shen, L. Hybrid hierarchical trajectory planning for a fixed-wing UCAV performing air-to-surface multi-target attack. J. Syst. Eng. Electron., 2012, 23(4), 256-264.[Full text via CrossRef]

9.     Besada-Portas, E.; Torre, L.d.l.; Cruz, J.M.d.l. & Andres-Toro, B.d. Evolutionary trajectory planner for multiple UAVs in realistic scenarios. IEEE Trans.. Robot., 2010, 26(4), 619-634.[Full text via CrossRef]

10.   Stevens, B.L. & Lewis, F.L. Aircraft control and simulation. John Wiley & Sons, Inc. Ed 2, New Jersey, 2003. 680 p.

11.   Lewis, L.-P.R. Rapid motion planning and autonomous obstacle avoidance for unmanned vehicles. Naval Postgraduate School, Monterey, California, USA, 2006. MA Thesis.

12.   Michel, F. Flatness and defect of nonlinear systems: Introductory theory and examples. Int. J. Control, 1995, 61(6), 1327-1361.[Full text via CrossRef]

13.   Zhang, Y.; Chen, J. & Shen, L. Real-time trajectory planning for UCAV air-to-surface attack using inverse dynamics optimization method and receding horizon control. Chinese J. Aeronaut., 2013, 26(4), 1038-1056.[Full text via CrossRef]

14.   Foo, J.L.; Knutzon, J.; Kalivarapu, V.; Oliver, J. & Winer, E. Path planning of unmanned aerial vehicles using B-splines and particle swarm optimization. J. Aeros. Comp. Inf. Com., 2009, 6, 271-290.[Full text via CrossRef]

15.   E. Rutquist, P. & M. Edvall, M. PROPT - Matlab optimal control software. TOMLABoptimization. 2010.

Dr Xueqiang Gu is a PhD scholar in College of Mechatronic Engineering and Automation, National University of Defense Technology (NUDT), Changsha, Hunan, China. His current research interests included: combat vehicle mission planning and intelligence control.

Dr Yu Zhang received his PhD in control science and engineering from NUDT, Changsha, Hunan, China in 2012. He is a Docent of College of Mechatronic Engineering and Automation in NUDT. His current research interests included: combat vehicle mission planning and artificial intelligence

Dr Jing Chen obtained his PhD in control science and engineering from NUDT in 1999. He is a Professor of College of Mechatronic Engineering and Automation in NUDT, Changsha, Hunan, China. His current research interests included: artificial intelligence and mission planning of aircraft

Mr Lincheng Shen is the Dean and Professor of College of Mechatronic Engineering and Automation in NUDT, Changsha, Hunan, China. He has published over 100 technical papers in refereed international journals and academic conferences. His current research interests included: mission planning, SAR image processing, biomimetic robotics, automation and control engineering.