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

    c = c + AMPLmatrix(fid,'obst',obst);

end

% defining of number of waypoints

nTarget = size(xT,1);

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

% defining of fuel penalty

FP=0.001/(nVeh*max(Amax)*T);

c = c + AMPLscalar(fid,'FP',FP);

% defining of individual completion time penalty

ITP=0.01/nVeh;

c = c + AMPLscalar(fid,'ITP',ITP);

% big M for logical relaxations

bigM=1000;

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

% defining of initial states

State0 = x0;

c = c + AMPLmatrix(fid,'State0',State0);

% defining of target points

StateF = xT';

c = c + AMPLmatrix(fid,'StateF',StateF);

% defining of circle approximation for speed limit

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

c = c + AMPLscalar(fid,'Crs',1);

% defining of speed limits

spdlimit=zeros(nVeh,2,Ncs);

for vv=[1:nVeh],

  Crs = Vmax(vv)*cos(0.5*alphas);

  for ii=[1:Ncs],

    spdlimit(vv,1,ii)=(1/Crs)*sin(ii*alphas);

    spdlimit(vv,2,ii)=(1/Crs)*cos(ii*alphas);

  end

end

c = c + AMPLmtrx3d(fid,'spdlimit',spdlimit,[1:nVeh],[1:2],[1:Ncs]);

% defining of circle approximation for acceleration limit

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

c = c + AMPLscalar(fid,'Crf',1);

% defining of force limits

frclimit = zeros(nVeh,2,Ncf);

for vv=[1:nVeh],

  Crf=Amax(vv)*cos(0.5*alphaf);

  for ii=[1:Ncf],

    frclimit(:,1,ii)=(1/Crf)*sin(ii*alphaf);

    frclimit(:,2,ii)=(1/Crf)*cos(ii*alphaf);

  end

end

c = c + AMPLmtrx3d(fid,'frclimit',frclimit,[1:nVeh],[1:2],[1:Ncf]);

% check for terminal region parameters

if isempty(targetBox)

  targetBox = zeros(nTarget,2);

end

c = c + AMPLmatrix(fid,'TargBox', targetBox');

% target region penalty

if isempty(targetPenalty)

  targetPenalty = zeros(nTarget,2);

end

c = c + AMPLmatrix(fid,'TargPen', targetPenalty');

% capability matrix

if isempty(capab)

  capab = ones(nVeh,nTarget);

end

c = c + AMPLmatrix(fid,'capab', capab);

% avoidance box

if isempty(vehBox)

  vehBox = zeros(nVeh,1);

end

c = c + AMPLvector(fid,'VehBox', vehBox);

% finished writing data file

fclose(fid);

% run AMPL

!ampl Trajectory.run

% load data files

load XOPT.dat

load FTOPT.dat

% plotting

close all

figure

hold on

% parameters

nVeh = length(FTOPT);

nStep = size(XOPT,1)/nVeh;

% plot obstacles

for ii = 1:size(obst,1),

   plot(obst(ii,[1 4 3 1 1]),obst(ii,[2 2 4 3 2]),'linewidth',5)

end

% plot targets

StateF=xT';

for ii=[1:size(StateF,2)],

    plot(StateF(1,ii),StateF(2,ii),'ro','markersize',24,'linewidth',2)

    text(StateF(1,ii),StateF(2,ii)-.3,['Final position']);

end

plotcols=['ro';'co';'bo';'go']; %colours on the plot

% shows only the segment before reaching the goal

for nn=[1:nVeh],

   tF=FTOPT(nn);

   pRng=[(nStep*(nn-1) + 1):(nStep*(nn-1) + 1 + tF)] ;

   plot(XOPT(pRng,1), XOPT(pRng,2), ...

        plotcols(nn,:),'markersize',3,'linewidth',6)

text(XOPT(pRng(1),1)-.2, XOPT(pRng(1),2)-.3,['Initial position']);

end

  axis equal %specifying of equal axis on the graph

3. Obtained graphical representations of work done

On the figures 4-2 to 4-12 the results of modeling are presented. During the modeling different initial data were used: speed of UAV, its turn rate, number of final positions, number of obstacles that must be overcome and various geometry of the obstacles. Unfortunately, because of usage of student version of AMPL CPLEX program module it was impossible to represent all possibilities of the program.

Figure 4-2: One obstacle and one final position.

Figure 4-3: One obstacle and two final positions.

Figure 4-4: One obstacle and two final positions.

Figure 4-5: Two obstacles and one final position.

Figure 4-6, 4-7: The same initial and final positions, different turn rate and speed

 

Figure 4-8, 4-9: The same initial and final positions but different trajectory because of different turn rates.

Figure 4-10: Final position between two obstacles.

 

Figure 4-11, 4-12: Obstacles of different geometry.

    

 

 

 

CONCLUSIONS

The selection of the goal is expressed in MILP form using binary variables that include other logical constraints such as collision avoidance. The point mass model captures the essential features of the aircraft dynamics such as speed and minimum turning radius constraints. The detailed analysis also clarified the effect of the time discretization. The non-convex minimum speed constraints have also been added to the problem to better capture the vehicle behavior in highly constrained environments. The effect on the computation time was also examined.

Experiments have been presented to demonstrate the use of Mixed-Integer Linear Programming for on-line replanning to control vehicles in the presence of dynamic uncertainty.

This initial set of experiments will be extended to involve more complicated scenarios. This work will include the distribution of the MILP planning process across multiple computers and transition to a greater number and variation of vehicles.

 

 

 

 

 

 

 

 

REFERENCES

1. Yoshiaki Kuwata. Real-time Trajectory Design for Unmanned Aerial Vehicles using Receding Horizon Control. – 2003. –Massachusets institute of technology.
2. Tom Schouwenaars, Bart De Moor. Mixed integer programming for multi-vehicle path planning // ECC2001 Conference. – Cambridge:2001
3. Arthur Richards, Yoshiaki Kuwata, Jonathan How. Experimental Demonstrations of Real-time MILP Control // AIAA Guidance, Navigation, and Control Conference and Exhibit 11-14 August 2003, Austin; AIAA 2003-5802
4. Tom Schouwenaars, Mario Valenti, Eric Feron, Jonathan How.Implementation and Flight Test Results of MILP-based UAV Guidance, Cambridge:2005
5. ILOG CPLEX 9.0 User’s guide. ILOG: 2003.
6. Ben Gardiner, Waseem Ahmad, Travis Cooper. Collision Avoidance Techniques for Unmanned Aerial Vehicles Technical Report #CSSE11 – 01, Auburn:2011
7. Yoshiaki Kuwata .Trajectory Planning for Unmanned Vehicles using Robust Receding Horizon Control
8. MATLAB, http://en.wikipedia.org/wiki/MATLAB.