Mixed-integer linear programming collision avoidance method in unmanned aerial vehicles trajectory planning

 

ABSTRACT

The explanatory note to the term paper “Mixed-integer linear programming collision avoidance method in unmanned aerial vehicles trajectory planning”:

The object of investigation – unmanned aerial vehicles.

The objective of the work – modeling of the collision avoidance scenario on the base of mixed-integer linear programming method.

Unmanned aerial vehicles (UAV) are widely used today. The investigation of optimal methods of collision avoidance will allow decreasing the probability of conflicts between UAVs and maintain the required safety level. The mixed-integer linear programming method is quite simple and useful for this purpose. It can be easily implemented not only in UAVs, but also in new systems of collision avoidance.

MIXED-INTEGER LINEAR PROGRAMMING, MILP, UNMANNED AERIAL VEHICLE, COLLISION AVOIDANCE, CONFLICT RESOLUTION, FLIGHTS SAFETY.

 

 

 

 

 

 

 

CONTENTS

LIST OF ACRONYMS…………………………………………………………..….5

INTRODUCTION………………………………………………………..………….6

1. PROBLEM FORMULATION………………………………………………….8

1. Trajectory Design………………………………………………..8 
2. Task Allocation…………………………………………………10

2. RECEDING HORIZON CONTROL USING MIXED-INTEGER LINEAR PROGRAMMING……………………………………………………………..11

1. Overview of the Receding Horizon Controller………………....11
2. Model of Aircraft Dynamics……………………………..……..12

1. Discrete Time System…………………………..…12
2. Speed Constraints…………….…………………..13 
3. Minimum Turning Radius………………….…….16 

3. Collision Avoidance Constraints………………………………..16

1. Obstacle Avoidance………….……………………17
2. Vehicle Avoidance……………………..…………18 

4. Plan beyond the Planning Horizon……………………………..19 
5. Target Visit..……………………………………………………21 
6. Cost Function……………………………………………………23

3. SOFTWARE……………………………………………………………………26
4. MODELING AND ITS RESULTS.……..……………………………………27

1. Annotation……………………………………………………...27
2. Text of the program…………………………………………….29
3. Obtained graphical representations of work done………………35

CONCLUSIONS………………………………………………………………..…..40

REFERENCES…………………………………………………………………..…41

 

LIST OF ACRONYMS

AMPL – A Mathematical Programming Language

CPLEX – Simplex method optimizer

MAC – Midair collision

MILP – Mixed-integer linear programming

MMKP  – Multiple-choice multi-dimension knapsack problem

MPC  – Model predictive control

NTG – Nonlinear trajectory generator

RHC – Receding horizon control

TOA – The time-step of arrival

UAV – Unmanned aerial vehicles

 

 

 

INTRODUCTION

Aviation safety is the number one priority for everyone working in the aviation industry. Over the past decade, the aviation community has witnessed a fundamental shift in its approach to safety. The evolution of the strategies is critical to ensure that international civil aviation remains the safest mode of transportation even as it continues to accommodate the significant growth in global populations and air travel forecast for the near future. Despite improvements in the general aviation safety record in recent years, the number of midair collisions (MACs) shows no corresponding decline. MACs continue to occur about 12 times a year on average, often resulting in multiple fatalities. That is why the lot of mathematical methods of collision avoidance was developed during the last decades and continues to improve day by day.

The object of this work is unmanned aerial vehicles (UAV) that is widely used all around the world. The role of UAV has changed over the past five to ten years as the level of autonomy onboard the vehicle increases. The primary advantage of UAVs is that they can perform missions in dangerous environments, such as battle fields and devastated areas, without risking the life of the human pilots. The variety of application areas include war zones, surveillance, forest fire monitoring, disaster relief, and weather forecast. Because of the absence of the main detector of probable conflict on UAV, it means human eyes, the strong and reliable method of conflict resolution must be used.

The objective of this work is to use mathematical programming techniques to find the optimal path between two states for a single vehicle. This path will be optimized with respect to both fuel and/or time, and must ensure that the vehicle does not collide with fixed obstacles. For this purpose the mixed-integer linear programming is used. It is well-suited for trajectory planning because it can directly incorporate logical constraints such as obstacle avoidance and waypoint selection and because it provides an optimization framework that can account for basic dynamic constraints such as turn limitations and maximum rate of climb.  Background The capabilities and roles of unmanned aerial vehicles (UAVs) are evolving, and new concepts are required for their control. For example, today’s UAVs typically require several operators per aircraft, but future UAVs will be designed to make tactical decisions autonomously and will be integrated into coordinated teams to achieve high-level goals, thereby allowing one operator to control a group of vehicles. Thus, new methods in planning and execution are required to coordinate the operation of a fleet of UAVs.

An overall control system architecture must also be developed that can perform optimal coordination of the vehicles, evaluate the overall system performance in real time, and quickly reconfigure to account for changes in the environment or the fleet. This thesis presents results on the guidance and control of fleets of cooperating UAVs, including goal assignment, trajectory optimization, and hardware experiments. For many vehicles, obstacles, and targets, fleet coordination is a very complicated optimization problem where computation time increases very rapidly with the problem size.  

 

 

 

 

 

 

 

1. PROBLEM FORMULATION
1. Trajectory Design 

Optimizing a kinematically and dynamically constrained path is a significant problem in controlling autonomous vehicles, and has received attention in the fields of robotics, undersea vehicles, and aerial vehicles. Planning trajectories that are both optimal and dynamically feasible is complicated by the fact that the space of possible control actions in extremely large and non-convex, and that simplifications reducing the dimensionality of the problem without losing feasibility and optimality are very difficult to achieve. Previous work demonstrated the use of mixed-integer linear programming (MILP) in off-line trajectory design for vehicles under various dynamic and kinematic constraints.

MILP allows the inclusion of non-convex constraints and discrete decisions in the trajectory optimization. Binary decision variables allow the choice of whether to pass “left” or “right” of an obstacle, for example, or the discrete assignment of vehicles to targets to be included in the planning problem. Optimal solutions can be obtained for these trajectory generation problems using commercially available software such as CPLEX. Using MILP, however, to design a whole trajectory with a planning horizon fixed at the goal in very difficult to perform in real time because the computational effort required grows rapidly with the length of the route and the number of obstacles to be avoided. This limitation can be avoided by using a receding planning horizon in which MILP in used to form a shorter plan that extends towards the goal, but does not necessarily reach it. This overall approach is known as either model predictive control (MPC) or receding horizon control (RHC) [1].

The performance of a RHC strongly depends on the proper evaluation of the terminal penalty on the shorter plan. This evaluation is difficult when the feasibility of the path beyond the plan must be ensured. Previous work presented a heuristic to approximate the trajectory beyond the shorter plan that used straight line paths 18 to estimate the cost-to-go from the plan’s end point to the goal. This RHC makes full use of the future states predicted through a model in order to obtain the current control inputs. In a trajectory design problem, the future states can be predicted by the straight line paths. This is because the shortest path from any location to the goal is a connection of straight lines, if no vehicle dynamics are involved and no disturbances act on the vehicle. As such, generating a coarse cost map based on straight line approximations (and then using Dijkstra’s algorithm) provides a good prediction of the future route beyond the planning horizon. The receding horizon controller designs a detailed trajectory over the planning horizon by evaluating the terminal state with this cost map. While this planner provides good results in practice, the trajectory design problem can become infeasible when the positioning of nearby obstacles leaves no dynamically feasible trajectory from the state that is to be reached by the vehicle. This is because the Dijkstra’s algorithm neglects the vehicle dynamics when constructing the cost map. In such a case, the vehicle cannot follow the heading discontinuities where the line segments intersect since the path associated with the cost-to-go estimate is not dynamically feasible. 

This work presents one way to overcome the issues with the previous RHC approach by evaluating the cost-to-go estimate along straight line paths that are known to be “near” kinodynamically feasible paths to the goal. This new trajectory designer is shown to: (a) be capable of designing trajectories in highly constrained environments; and (b) travel more aggressively since the turn around an obstacle corner is designed at the MILP optimization phase, unlike a trajectory designer that allows only predetermined turns at each corner.

Although RHC has been successfully applied to chemical process control, applications to systems with faster dynamics, such as those found in the field of aerospace, have been impractical until recently. The on-going improvements in computation speed has stimulated recent hardware work in this field, including the on-line use of nonlinear trajectory generator (NTG) optimization for control of an aerodynamic system. The use of MILP together with RHC enables the application of the real-time trajectory optimization to scenarios with multiple ground and aerial vehicles that operate in dynamic environments with obstacles.  

2. Task Allocation

 

 A further key aspect of the UAV problem is the allocation of different tasks to UAVs with different capabilities. This is essentially a multiple-choice multi-dimension knapsack problem (MMKP), where the number of possible allocations grows rapidly as the problem size increases. The situation is further complicated if the tasks: - are strongly coupled – e.g., a waypoint must be visited three times; -have tight relative timing constraints – e.g., three UAVs must be assigned to strike a target from three different directions within 2 seconds of each other [2]. 

With limited resources within the team, the coupling and timing constraints can result in very tight linkages between the activities of the various vehicles. Especially towards the end of missions, these tend to cause significant problems (e.g., “churning” and/or infeasible solutions) for the approximate assignment algorithms based on “myopic algorithms” (e.g. iterative “greedy” or network flow solutions) that have recently been developed. MILP, again, provides a natural language for codifying these various mission objectives and constraints using a combination of binary (e.g., as switches for the discrete/logical decisions) and continuous (e.g., for start and arrival time) variables. MILP allows us to account for the coupling and relative timing constraints and to handle large fleets with many targets and/or pop-up threats. 

 

 

2. RECEDING HORIZON CONTROL USING MIXED-INTEGER LINEAR PROGRAMMING

2.1 Overview of the Receding Horizon Controller 

Improvements in UAV capabilities make it possible for UAVs to perform longer and more complicated missions and scenarios. In these missions and scenarios, optimal path navigation through complicated environments is crucial to mission success. As more vehicles and more targets are involved in the mission, the complexity of the trajectory design problem grows rapidly, increasing the computation time to obtain the optimal solution. One alternative to overcome this computational burden is to use receding horizon control (RHC).

The RHC uses a plant model and an optimization technique to design an input trajectory that optimizes the plant’s output over a period of time called the planning horizon. Ap ortion of the input trajectory is then implemented over the shorter execution horizon, and the optimization is performed again starting from the state that is to be reached. If the control problem is not completed at the end of the planning horizon, the cost incurred past the planning horizon must be accounted for in the cost function. The selection of the terminal penalty in RHC design is a crucial factor in obtaining reasonable performance, especially in the presence of obstacles and no-fly zones [1].  

In general, the cost function of a receding horizon controller’s optimization problem estimates the cost-to-go from a selected terminal state to the goal. For vehicle trajectory planning problems in a field with no-fly zones, presented a receding horizon controller that uses the length of a path to the goal made up of straight line segments as its cost-to-go. This is a good approximation for minimum time of arrival problems since the true minimum distance path to the goal will typically touch the corners of obstacles that block the vehicle’s path. In order to connect the detailed trajectory designed over the planning horizon and the coarse cost map beyond it, the RHC selects an obstacle corner that is visible from the terminal point and is associated with the best path. This approach has another advantage in terms of real-time applications.

The cost map not only gives a good prediction of vehicle behavior beyond the planning horizon, but because it is very coarse, in also can be rapidly updated when the environment and/or situational awareness changes. The following sections focus on solving the vehicle trajectory design problem after a set of ordered goals are assigned to each vehicle. The MILP formulation presented here extends an existing algorithm to incorporate more sophisticated scenarios (e.g., multiple vehicles, multiple goals) and detailed dynamics (e.g., constant speed operation). 

2.2 Model of Aircraft Dynamics 

It has been shown that the point mass dynamics subject to two-norm constraints form in approximate model for limited turn-rate vehicles, provided that the optimization favors the minimum time, or minimum distance, path. By explicitly including minimum speed constraints, this formulation can be applied to problems with various objectives, such as expected score and risk, with a minimal increase in computation time. 

2.2.1. Discrete Time System 

Aircraft dynamics is expressed as a simple point mass with position and velocity     [x, y, v_x, v_y]^T as state variables and acceleration [a_x, a_y]^T as control inputs [1],[2].

(2.1)

The zero-order hold equivalent discrete time system in

(2.2)

     where k expresses a time step and Δt is the time interval. Note that the control input [a_x, a_y]_k ^T stays constant over each time interval Δt under the zero-order hold assumption. 

2.2.2 Speed Constraints 

The constraint on the maximum speed v_max is written as a combination of linear constraints on the velocity vector v using uniformly distributed unit vectors

	(2.3)
	(2.4)

where n_vmax is the number of unit vectors i_k onto which the velocity vector in projected. The speed constraint is effectively a constraint on the length of the projection of the velocity vector onto a unit vector. Eqs. 2.3 and 2.4 require that the velocity vector be inside a regular polygon with n_vmax sides circumscribed about a circle of radius v_max. Aconstrain t on the minimum speed v_min can be expressed in the similar way: the dot product of the velocity vector and a unit vector must be larger than v_min. However, it in different from the maximum speed constraint in that at least one of the constraints must be active, instead of all of them,  

	(2.5)
	(2.6)

        Here, n_vmin is the number of unit vectors onto which the velocity vector is projected. Eq. 2.5 is a non-convex constraint and is written as a combination of mixed-integer linear constraints

(2.7)

 

(2.8)

        where M_v is a number larger than 2v_min, and the b_speed,k are binary variables that express the non-convex constraints in the MILP form.

Note that if b_speed,k = 1, the inequality in Eq. 2.7 is active, indicating that the minimum speed constraint is satisfied. On the other hand, if b_speed,k = 0, the k^th constraint in Eq. 2.7 is not active, and the constraint on minimum speed is relaxed. Eq. 2.8 requires that at least one constraint in Eq. 2.7 be active. Eqs. 2.6 to 2.8 ensure that the tip of the velocity vector lies outside of a regular polygon with n_vmin sides circumscribed on the outside of a circle with radius v_min.  Since the minimum speed constraints require binary variables, which typically slows down the MILP optimization process, n_vmin is typically much smaller than n_vmax . Note that the speed constraints can become infeasible if v_max and v_min are set equal unless n_vmax = n_vmin.   

2.2.3 Minimum Turning Radius 

UAVs usually fly at roughly a constant speed v and have a minimum turning radius r_min. The constraint on the turning radius r is 

(2.9)

         and this constraint can be written as a constraint on lateral acceleration,

(2.10)

        where a is the magnitude of the acceleration vector and a_max is the maximum acceleration magnitude. Similar to the maximum speed constraint, the constraint on the maximum acceleration is written as a combination of linear constraints on the acceleration vector a 

(2.11)

 

(2.12)

       where n_a is the number of unit vectors onto which the acceleration vector is projected. The constraints on velocity in Eqs. 2.3 to 2.8 keep the speed roughly constant, so the allowable acceleration vector direction is perpendicular to the velocity vector [3]. The state equation (Eq. 2.2) for the vehicle uses a zero-order hold on the inputs, so the acceleration vector stays the same over the time-step Δt. Figure 2-3 shows a turn with the maximum lateral acceleration. The actual minimum turning radius r_{min actual} is a radius of the polygon. It is geometrically calculated, using the relation between a_max and r_min in Eq. 2.10, as

(2.13)

       and is slightly smaller than r_min.  

 

2.3. Collision Avoidance Constraints 

2.3.1. Obstacle Avoidance 

During the overall mission, UAVs must stay outside of no-fly-zones, which are modelled as obstacles in our formulation. Obstacle avoidance is a non-convex constraint and requires binary variables and a large number M in MILP. At each time step k, the vehicle must stay outside of a rectangular obstacle defined by four parameters [x_l, y_l, x_u, y_u].

The four edges of each obstacle are given by 

x = xl; x = xu; y = yl; or y = yu

(2.14)

The constraints on vehicle position are formulated as 

xk ≤ xl +M bobst,1jk yk ≤ yl +M bobst,2jk xk ≥ xu −M bobst,3jk yk ≥ yu −M bobst,4jk   j=1, …, no, k=1, …, np

(2.15) (2.16) (2.17) (2.18) (2.19)

       where [x_k, y_k]^T gives the position of a vehicle at time k, M is a number larger than the size of the world map, n_o is the number of obstacles, n_p is the number of steps in the planning horizon, and b_obst,ijk is an i − j − k^th element of a binary matrix of size 4 by n_o by n_p. At time-step k, if b_obst,ijk = 0 in Eqs. 2.15 to 2.18, then the constraint is active and the vehicle is outside of the j^th obstacle. If not, the large number M relaxes the obstacle avoidance constraint.

The logical constraint in Eq. 2.19 requires that at least one of the four constraints be active for each obstacle at each time-step. The obstacle shape is assumed to be rectangular, but this formulation is extendable to obstacles with polygonal shapes. Also, non-convex obstacles can be easily formed by overlapping several rectangular obstacles. Figure 2-1 shows that we must expand the obstacle size at the planning level to account for the discrete steps taken by the vehicle. This increase is done at both the estimationand trajectory design phases. Since the avoidance constraints are only applied at discrete time steps shown as ⊗ marks, it is possible for the planned trajectory to “cut the corner” of the obstacle between time points. Each waypoint is separated by vΔt and an obstacle [x_l, y_l, x_u, y_u] must be expanded in each direction by , which is the maximum incursion distance, so that

(2.20)

Figure 2-1: Corner cut and obstacle expansion

2.3.2. Vehicle Avoidance 

Collision avoidance between vehicles is written in the same way as obstacle avoidance. Assume the i^th vehicle has a certain physical size and safety distance surrounding in that forms together a rectangle of 2d_i by 2d_i around its center. At each time k, the position of the i^th vehicle and that of the jt^h vehicle must satisfy the following relations [3]: 

xik ≤ xjk + (di + dj) +M bveh,1ijk yik ≤ yjk + (di + dj) +M bveh,2ijk   xik ≥ xjk − (di + dj) −M bveh,3ijk yik ≥ yjk − (di + dj) −M bveh,4ijk   i = 1, . . . , nv,   j= i + 1, . . . ,  nv     k = 1, . . . , np

(2.21) (2.22) (2.23) (2.24) (2.25)

where [x_ik, y_ik]^T gives the position of a vehicle i at time k, n_v is the number of vehicles, and b_veh,_lijk is the l − i − j − k^th element of a binary matrix of size 4 by n_v by n_v by n_p. If b_veh,_lijk = 0 in Eqs. 2.21 to 2.24, the constraint is active, and the safety boxes of the two vehicles, i and j, do not intersect each other. Again, in order to account for the discretized time, a margin is added to the vehicle rectangle of 2d_i by 2d_i. If two vehicles i and j are moving towards each other with a speed of v_i and v_j respectively, the distance between each waypoint in a relative coordinate frame, can be as large as (v_i + v_j)Δt. Thus, the size of the ith vehicle in this MILP optimization model must be expanded to a size of 2d_i|expanded by 2d_i|expanded where

(2.26)

2.4. Plan beyond the Planning Horizon 

The Receding Horizon Controller uses a coarse cost map based on straight lines to predict the trajectory beyond the planning horizon. There are two aspects involved in connecting the detailed trajectory over the planning horizon to this coarse cost map. First, MILP selects a cost point that leads the vehicle to the goal along the shortest path. Second, it ensures that the selected cost point is visible from a point on the planning horizon (i.e., the straight line segment connecting the selected cost point and the horizon point must be obstacle/collision free). Given the location of the obstacles and a goal, the cost points are defined as the obstacle corners and the goal itself. The shortest distance from a cost point to the goal along kinematically feasible straight-line segments forms a cost associated with the cost point and goal. If a cost point is inside of another obstacle, it has an infinite cost. Let [x_cp,i, y_cp,i]^T denote the it^h cost point, c_i the cost associated with the i^th cost point, and c_vis,k the cost-to-go at the cost point selected by vehicle k. The binary variables associated with cost point selection b_cp will have three dimensions for cost point, goal, and vehicle, where n_c is a number of cost points, and n_g is a number of goals [4]:

(2.27)   (2.28)

Eq. 2.28 enforces the constraint that each vehicle must choose a combination of goal and cost point, and Eq. 2.27 extracts the cost-to-go at the selected point from the cost map c_i. In order to ensure that the selected cost point [x_vis,k, yv_is,k]^T is visible from the terminal point [(x_np)k, (y_np)k]^T of vehicle k, obstacle avoidance is checked at n_t test points that are placed on a line segment connecting the two points. This requires binary variables b_v  is that have four dimensions: the obstacle corner, vehicle, test point, and obstacle. The test conditions can be written as:

	(2.29)   (2.30)         (2.31)
k=1, …, nv, m=1, …, nt , n=1, …, no	(2.32) (2.33) (2.34) (2.35) (2.36)

where [x_LOS,k, y_LOS,k]^T in Eq. 2.30 is a line-of-sight vector from the terminal point to the selected cost point for vehicle k, and the binary variable b_vis has four dimensions (obstacle boundary, vehicle, test point, and obstacle). 

2.5. Target Visit 

The highest goal of the mission is to search, attack, and assess specific targets. These tasks can be generalized as visiting goals and performing the appropriate action. The heading direction at the target can be included if assessment from a different angle is required. In this section, these goals are assumed to be allocated among several vehicles so that each vehicle has an ordered list of goals to visit. The goals are assumed to be set further apart than the planning horizon, so each vehicle can visit at most one goal point in each plan. The constraint on the target visit is active at only one time-step if a plan reaches the goal, and is relaxed at all of the waypoints if it does not. If a vehicle visits a goal, the rest of the steps in the plan are directed to the next goal.

Therefore, only the first two unvisited goals are considered in MILP (i.e., the number of goals n_g is set to be 2). The selection of the cost point and the decision of the time of arrival are coupled as 

(2.37)   (2.38)

where b_arrival,ij = 1 if the j^th vehicle visits the first target at the i^th time step, and equals 0 if not. Also, i = n_t + 1 indicates the first target was not reached within the planning horizon. If a vehicle cannot reach the goal or reaches the goal at the final step, Eq.2.37 must use the cost map built for the first goal. If a vehicle can visit the first goal in the planning horizon, Eq. 2.38 requires it to use the cost map built about the second goal. Agoal has a v_iΔt by v_iΔt square region of tolerance around it because the system in operated in discrete time. This tolerance is encoded using the intermediate variable

derr          i=1, …, nt, j = 1, …, nv .

(2.39) (2.40) (2.41) (2.42)   (2.43)

If a vehicle j reaches its goal at time-step i, then b_arrival,ij = 1 and the constraints in Eqs. 2.39 to 2.42 are active; if not, all the constraints are relaxed. Including d_err in the objective function causes the point of visit x_ij (where b_arrival,ij = 1) to move closer to the goal. The heading constraints can be written in the same way 

2.6. Cost Function 

This formulation minimizes the time of arrival as a primary objective. Under the assumption that each vehicle flies at a constant speed, the minimization of distance is equivalent to the minimization of time. The distance to the goal is made up of three segments: 1) initial position to the terminal point of the plan; 2) terminal point to the cost point; and 3) the precalculated cost-to-go at the selected cost point. Since the distance from the initial position to the terminal point is constant, the MILP optimizer minimizes the sum of the length of the line-of-sight vector 2) and the cost-to-go at the selected cost point 3) by choosing the best combination of terminal and cost-to-go points [3].

The length of the line-of-sight vector is a two norm of the line-of-sight vector defined in Eq. 2.30. This is calculated in the MILP by minimizing an upper bound:

(2.44)   (2.45)

where n_l is the number of unit vectors onto which the line-of-sight vector is projected and l_j gives an upper bound of the line-of-sight vector length for vehicle j. In this case, the objective function J₀ to be minimized in

(2.46)

 

 where c_vis,j is a cost-to-go at the cost point selected by the jth vehicle as defined in Eq. 2.28. The time-step of arrival (TOA) for the jth vehicle is expressed as

(2.47)

Note that in Eq. 2.47, if k = n_t+1 when b_arrival,kj = 1, then the j^th vehicle will not reach the target in the planning horizon, thereby resulting in a higher cost. Combined with Eq. 2.46, the term TOA forms the objective function for the minimum time of arrival problem:

(2.48)

where v_jΔt converts the discrete time-step (TOA_j is an integer) into distance and α is a large weighting factor that makes the first term dominant. With a small α, the vehicle tends to minimize the cost-to-go, so that it arrives at the first goal at the terminal step in the planning horizon. This is in contrast to the more desirable case of arriving at the first goal before the terminal step, and then minimizing the cost-to-go to the next goal. Thus the weight α must be larger than the number of steps required to go from the first goal to the next goal.  Finally, the term J_a, containing two auxiliary terms, is added to the primary cost in Eq. 2.48 to decrease the position and heading error at the target. The new cost to be minimized is defined as

(2.49)

where

(2.50)

Note that an overly large value of β results in “wobbling” near the target. The wobbling occurs as the vehicle attempts to maneuver, such that on the same time-step when it enters the target region, it arrives at the center of the target. Also note that the equations involving the dynamics, constraints, and costs are not dependent of the number of vehicles except for vehicle collision avoidance, and thus these can be solved separately if vehicle collision avoidance is not an issue [7]. 

 

 

 

 

 

 

 

 

 

3. SOFTWARE

The total optimization problem from the previous section can be readily solved by a mixed integer program solver implemented in the CPLEX software package [5]. CPLEX requires the standard MPS-format as input, which can be generated by running a specific m-script in MATLAB. However, writing such an MPS-file generating script is an error prone task and such scripts are usually not very flexible in adapting the cost function or the type of constraints. I have therefore used the A Mathematical Programming Language (AMPL) to formulate the path planning problem. Implementing the constraints in AMPL is straightforward, requiring minimum translation. The forms of the problem and the constraints are defined in a model file, while the parameter values are in a separate data file.

Therefore, changes to the problem may be made without rebuilding the constraints expressions. The data files can be edited directly or generated by a simple MATLAB script. AMPL combines the model and data to an MPS-file, which is then solved by CPLEX. Moreover, AMPL can handle piece-wise linear variables and thus absolute values: it hence generates the necessary slack variables and extra constraints automatically. I thus specify the system matrices, the obstacle coordinates and other optional dynamical parameters in MATLAB, and use a specific script to generate an AMPL input file. An AMPL script selects the appropriate model and data file, calls the solver and writes the resulting solution to an output file readable by MATLAB. Other MATLAB -scripts then plot the path and visualize the state and input sequence.    

 

 

 

4. MODELING AND ITS RESULTS
1. Annotation

Introductory passage. In accordance to previously decribed method the program “Modeling of unmanned aerial vehicle trajectory on the base of MILP” was created in programming environment MATLAB version 7.10.0.499. As it was mentioned, in order to implement MILP, the additional complex of computer program was used, namely AMPL CPLEX [5] software optimization package, that is an algebraic modeling language for describing and solving high-complexity problems for large-scale mathematical computation.

CPLEX solves linear and convex quadratic programs by simplex or interior-point methods, and linear and convex quadratic integer programs by a branch-and-bound procedure.

Functional area. The main objective of the program is simulation of UAV trajectory planning and finding of optimal way from the initial to final point with a glance to time and fuel consumption with the overcoming of obstacles. The program can be used on all computers the operating systems of which are compatible with MATLAB. The implementation of the program requires the using of AMPL modules. During the creation of program the student version of AMPL (CPLEX) was used. That’s why there is a limitation on input data because the student edition is limited to 300 variables and 300 constraints.

Description of logic.The program consists of the following structural parts: initializing of structure, dynamics model specifying, environmental conditions specifying, solving of the problem by MILP, plotting of the input and output parameters.

 The program starts with the initialization of structure. Each field in the structure appears to be a separate MATLAB array. After that the parameters of dynamic model (speed, turn rate, time step and number of steps) and environmental model (initial state of UAV, its final position (one or more) and obstacle (one or more) must be inserted.

The next step after initial data entry is converting the described dynamic model into acceleration limits and definition of approximations for magnitude limits. After that the data file can be started to be written. AMPL allows defining the initial data and rewriting it into the new file. All previously specified and all following data is written to the files Trajectory.dat and Trajectory.mod. During the writing of these files the following steps of problem solving are implemented with the help of AMPL module: convertion of the system into 2-dimensional in AMPL format, defining of the number of timesteps, vehicles and obstacles, initial and final states which are specified at the beginning of the program. It is compulsory to define speed and force limits and circle approximation for acceleration limit. Also the significant part is the calculation of fuel and time completion penalty, which allow choosing the most fuel and time-effective way for UAV. At the end of the program after all required computations the plots of the obstacles and UAV trajectories are created.

MATLAB (matrix laboratory) is a numerical computing environment and fourth-generation programming language. Developed by MathWorks, MATLAB allows matrix manipulations, plotting of functions and data, implementation of algorithms, creation of user interfaces, and interfacing with programs written in other languages [8].

The algorithm of program “Modeling of unmanned aerial vehicle trajectory on the base of MILP” is presented on the Fig 4-1:

 

 

 

 

 

 

 

 

 

 

 

 

 

Figure 4-1. The algorithm of mixed-integer linear programming model of UAV trajectory

2. Text of the program

The text of the program “Modeling of unmanned aerial vehicle trajectory on the base of MILP” on the initial language:

clear all %cleaning of command window and graphical windows

clc

clf

% initializing of structure

Nt             = [];

dt             = [];

Vmax           = [];

Wmax           = [];

x0             = [];

xT             = [];

obst           = [];

vehBox         = [];

targetBox      = [];

targetPenalty  = [];

capab          = [];

% specifying of dynamics models

Vmax=0.2;                % speed [m/s]

Wmax=(pi/180)*4;         % turn rate [rad/s]

dt=2;                    % time step [sec]

Nt=10;                   % number of steps

% specifying of environment

theta=50/180*pi;

x0=[-4.5 -3, Vmax*.5*[cos(theta) sin(theta)]];   % initial states

xT=[-1.5 -3; -3 -3; -2 -1];                      % target positions

obst=[ -4 -3 -2 -3];                             % obstacles

%----------------------------------------------------------------------

% solve MILP

%----------------------------------------------------------------------

% convert dynamics to acceleration limits

Amax=(Vmax).*(Wmax);

% approximations for magnitude limits

Ncf=8;

Ncs=8;

alphas=2*pi/Ncs;

alphaf=2*pi/Ncf;

% start of writing file

c=0;

fid=fopen('Trajectory.dat','w');

c = c + AMPLcomment(fid,'Matlab generated AMPL data file\n');

c = c + AMPLcomment(fid,'');

c = c + AMPLcomment(fid,'For use with model Trajectory.mod');

c = c + AMPLcomment(fid,'');

% convertion of all in two dimensions

c = c + AMPLscalarint(fid,'nDim', 2);

% defining of number of vehicles

nVeh = size(x0,1);

c = c + AMPLscalarint(fid,'N',nVeh);

% defining of number of timesteps

T=Nt;

c = c + AMPLscalarint(fid,'T',T);

% discrete system

Ad = [eye(2),   dt*eye(2);

      zeros(2),    eye(2)];

Bd = [0.5*dt^2*eye(2);

            dt*eye(2)];

% writing of system matrices

c = c + AMPLmatrix(fid,'A',Ad);

c = c + AMPLmatrix(fid,'B',Bd);

% size

[nState,nInput] = size(Bd);

c = c + AMPLscalarint(fid,'nState',nState);

c = c + AMPLscalarint(fid,'nInput',nInput);

% defining of number of obstacles

nObst = size(obst,1);

c = c + AMPLscalarint(fid,'nObst', nObst );

if nObst

    % obstacles