2023年深圳杯数学建模C题无人机协同避障航迹规划解题全过程文档及程序

2023年深圳杯数学建模

C题 无人机协同避障航迹规划

原题再现：

  平面上A、B两个无人机站分别位于半径为500 m的障碍圆两边直径的延长线上，A站距离圆心1 km，B站距离圆心3.5 km。两架无人机分别从A、B两站同时出发，以恒定速率10 m/s飞向B站和A站执行任务。飞行过程中两架无人机必须避开障碍圆、并且不得碰面（即两架无人机的连线必须保持与障碍圆处于相交状态）。无人机的转弯半径不小于30 m。请建立数学模型，解决以下问题：

  问题1 要求两架无人机中第一个到达目的站点的用时最少，给出两架无人机的飞行航迹方案。
  问题2 要求两架无人机中第二个到达目的站点的用时最少，给出两架无人机的飞行航迹方案。
  问题3 当B站点到圆心的距离变化(其他参数保持不变)时，问题1和问题2中的最优航迹会发生什么变化？
  问题4 当B机的恒定速率在[10,30] m/s内变化(其他参数保持不变)时，问题1和问题2中的最优航迹会如何变化？
  问题5 当B机的恒定速率在[10,50] m/s内变化、B站点到圆心的距离在[1,10] km内变化(其他参数保持不变)时，问题2中的最优航迹会如何变化？

整体求解过程概述(摘要)

  本文首先针对无人机必须跨过障碍圆且两架无人机的连线必须相交障碍圆，因此本文建立0-1多目标整数规划来解决无人机最优路径问题。并利用Dijkstra算法和蚁群算法来对目标模型进行求解，从而给出最佳航迹规划的方案。
  为了解决第一架无人机到达目的用时最少的航迹规划路径，本文建立0-1多目标整数规划模型，以无人机飞行路径最小和避开障碍圆为目标，从三个角度来考虑A、B站的位置，因其障碍圆位置并不固定，所以本文先求出点到点之间的距离，然后再结合障碍圆的约束条件，当A、B同向且在同一半径上，最短无人机路径长为（1.5+1000π）km，当A、B异向且在同一直径上，其最短路径为（3.5+1000π）km，当B在任意位置，其最短路径长为（1km/cosα+1000π-1）km。
  为了解决第二架无人机到达目的用时最少的航迹规划路径，本文根据第一问的分类讨论,采用多角度优化对无人机进行二维航迹的最短路径规划，即可以得到最佳路线。但是注意到A、B站的位置依旧任意，所以本文依旧采用分类讨论来对节点的选择上进行优化处理。因此，本文建模的整体思想是将上述无人飞机最短路径规划问题，本文采取多角度优化来对无人机航行位置和角度进行优化处理,当A、B同向且在同一半径上，最短无人机路径长为30.052arcsin2.376+30.002arcsin0.499km，当A、B异向且在同一直径上，其最短路径为30.169arcsin3.850+30.002arcs3.850km，当B在任意位置，其最短路径长为(〖S_min〗^2+3600)/120 arcsin (120*S_min)/(〖S_min〗^2+120)+30.002arcsin0.499km。为了解决当B站点到圆心的距离变化，第一架和第二架无人机最优路径的变化，本文采用了Dijkstra算法将无人机飞行的最短轨迹求解出来，并且在构建模型时采用了圆盘的方法。无人机在飞行过程中所产生的转弯半径强度大小就决定了无人机飞行过程中的转弯半径距离，在无人机飞行过程中所产生的转弯半径允许范围内,当无人机飞行时所转弯半径越近,所受到的威胁就越大。
  为了解决当B机无人机飞行的恒定速率控制在区间范围内对第一架和第二架无人机飞行的最优路径进行调整，因其本题中将恒定速率变成一个变量进行处理，那么必然导致其角度发生一定程度的变化，本文在第三问的基础上建立的转弯半径代价和转弯角度代价的基础上引入速率代价这个限制因素，因此本文选择蚁群算法来对两架无人机路径进行优化。当增加考虑无人机飞行过程中速率的增加和无人机飞行过程中距离的变化，第二架无人机最优航线的方案，本文将从两方面考虑，一方面在问题二中无人机在飞行过程中会有一个飞行角度的调整，另一方面随着B站点位置的随机性和无人机起飞速度的随机性，考虑飞行器在实际飞行过程中的30米的最小转弯半径，实现无人机航行轨迹的重新规划。使得航迹路程尽可能小的，同时满足其他限制条件，所以本文采用多目标动态规划解决问题。

模型假设：

  （1）忽略具体地势，天气和无人机操作等因素的影响。
  （2）不考虑无人机起飞和降落时的限制。
  （3）不考虑无人机燃料耗尽等因素对航迹的影响。
  （4）所规划的无人机航迹系统能够准确无误的执行。

问题分析：

  问题一的分析
  问题一要求两架无人机中第一个到达目的站点的用时最少，给出两架无人机的飞行航迹方案,因其题目中的障碍圆位置是随机的，而且无人机要求是不能碰面且两架无人机的连线要相交于障碍圆，从而本文利用二维平面，从而建立了多目标协同规划的航迹最短路径模型。因此要从三个方位因素来考虑A、B站的位置（这里主要指0°~180°），那么首先考虑A、B同向且在同一半径上，其次考虑A、B异向且在同一直径上，最后考虑B在任意位置，因其障碍圆位置并不固定，所以本文先求出点到点之间的距离，然后再结合障碍圆的约束条件，加上一个障碍圆的周长在减去直径的距离得出一个范围区间。
  问题二的分析
  问题二要求两架无人机中第二个到达目的站点的用时最少，本文根据第一问的分类讨论,采用多角度优化对无人机进行二维航迹的最短路径规划，可以得到最佳路线。但是注意到A、B站的位置依旧任意，所以本文依旧采用分类讨论来对节点的选择上进行优化处理。因此，综上所述本文建模的整体思想是将上述无人飞机最短路径规划问题，利用多角度优化的方法,对相应的无人机航行节点的选择进行优化处理,以便能在最短的时间内规划出无人机最优二维航迹路径。
  问题三的分析
  问题三要求改变B站到圆心的距离从而调整第一架和第二架无人机的最优路线，在无人机飞行过程中包含其无人机在航行过程中的转弯角度代价和无人机航行过程中的转弯半径代价。本文假设在相同转弯半径的作用下，无人机飞行航线也大致相同，因此，转弯半径强度大小就决定了转弯半径距离，其呈现为四次方的反比(1/d4)，本文认为,在转弯半径允许范围内,转弯半径越近,所受到的威胁就越大，反之则相反。因此本文借助圆盘的方法，构建出模型,并利用Dijkstra算法求解无人机飞行最优路径。
  问题四分析
  问题四要求将无人机飞行的恒定速率控制在一个区间范围内对第一架和第二架无人机飞行的路径进行调整优化，因其本题中将恒定速率变成一个变量进行处理，那么必然导致其角度发生一定程度的变化，在其他因素不变的情况下，本题自然而然的就变成了一个组合优化路径问题。由于同时考虑无人机的转弯半径、转弯速率和转弯角度的影响，本文在第三问的基础上建立的转弯半径代价和转弯角度代价的基础上引入速率代价这个限制因素，因此本文选择蚁群算法来对第一架和第二架无人机路径进行优化处理。
  问题五的分析
  问题五是在问题二的基础上，增加考虑无人机飞行过程中速率的增加，无人机飞行过程中距离的变化，从而导致其无人机飞行转弯角度约束。一方面在问题二中无人机在飞行过程中会有一个飞行角度的调整，另一方面随着题目B站点位置的随机性和无人机起飞速度的随机性，其必然会导致无人机飞行时的航迹变化，
  考虑飞行器在实际飞行过程中的30米的最小转弯半径，实现无人机航行轨迹的重新规划。使得航迹路程尽可能小的，同时满足其他限制条件，即为多目标动态规划问题。

模型的建立与求解整体论文缩略图

全部论文请见下方“ 只会建模 QQ名片” 点击QQ名片即可

程序代码：(代码和文档not free)

function [d,DD]=dijkstra(D,s)
%Dijkstra最短路算法Matlab程序用于求从起始点s到其它各点的最短路
%D为赋权邻接矩阵
%d为s到其它各点最短路径的长度
%DD记载了最短路径生成树
[m,n]=size(D);
d=inf.*ones(1,m);
d(1,s)=0;
dd=zeros(1,m);
dd(1,s)=1;
y=s;
DD=zeros(m,m);
DD(y,y)=1;
counter=1;
while length(find(dd==1))<m
     for i=1:m
         if dd(i)==0
             d(i)=min(d(i),d(y)+D(y,i));
         end
     end
     ddd=inf;
     for i=1:m
         if dd(i)==0&&d(i)<ddd
             ddd=d(i);
         end
     end
     yy=find(d==ddd);
     counter=counter+1;
     DD(y,yy(1,1))=counter;
     DD(yy(1,1),y)=counter;
     y=yy(1,1);
     dd(1,y)=1;
end

import numpy as np

import random

import math

import matplotlib.pyplot as plt

def aco_mtsp(points, atms_number, service_centers=2):

def fitness_pop(population):

fitness_result = np.zeros(len(population))

for i in range(len(fitness_result)):

fitness_result[i] = fitness(population[i])

return fitness_result

def fitness(creature):

sum_dist = np.zeros(len(creature))

for j in range(len(creature)):

mat_path = np.zeros((dist.shape[0], dist.shape[1]))

path = creature[j]

if len(path) != 0:

for v in range(len(path)):

if v == 0:

mat_path[engineers[j], path[v]] = 1

else:

mat_path[path[v - 1] + service_centers, path[v]] = 1 mat_path = mat_path * dist

sum_dist[j] = (np.sum(mat_path) + dist[engineers[j], path[-1]]) / velocity + repair_time * len(path)

return np.max(sum_dist)

def birth_prob(fitness_result):

birth_prob = np.abs(fitness_result - np.max(fitness_result))

birth_prob = birth_prob / np.sum(birth_prob)

return birth_prob

def mutate(creat, engi):

pnt_1 = random.randint(0, len(creat) - 1)

pnt_2 = random.randint(0, len(creat) - 1)

if random.random() < mut_1_prob:

creat[pnt_1], creat[pnt_2] = creat[pnt_2], creat[pnt_1]
if random.random() < mut_2_prob and pnt_1 != pnt_2:

if pnt_1 > pnt_2:

pnt_1, pnt_2 = pnt_2, pnt_1

creat[pnt_1:pnt_2 + 1] = list(reversed(creat[pnt_1:pnt_2 + 1]))

if random.random() < mut_3_prob:

engi = [number - 1 for number in engi if number != 0]

#	engi = [number - 2 for number in engi if number > 1] while (sum(engi) != atms_number):

engi[random.randint(0, len(engi) - 1)] += 1 return creat, engi

def two_opt(creature):

sum_dist = np.zeros(len(creature))

for j in range(len(creature)):

mat_path = np.zeros((dist.shape[0], dist.shape[1]))

path = creature[j]

if len(path) != 0:

for v in range(len(path)):

if v == 0:

mat_path[engineers[j], path[v]] = 1

else:


mat_path[path[v - 1] + service_centers, path[v]] = 1 mat_path = mat_path * dist

sum_dist[j] = (np.sum(mat_path) + dist[engineers[j], path[-1]]) / velocity + repair_time * len(path)

for u in range(len(creature)):

best_path = creature[u].copy()

while True:

previous_best_path = best_path.copy()

for x in range(len(creature[u]) - 1):

for y in range(x + 1, len(creature[u])):

path = best_path.copy()

if len(path) != 0:

path = path[:x] + list(reversed(path[x:y])) + path[y:] # 2-opt swap mat_path = np.zeros((dist.shape[0], dist.shape[1]))

for v in range(len(path)):

if v == 0:

mat_path[engineers[u], path[v]] = 1

else:

mat_path[path[v - 1] + service_centers, path[v]] = 1 mat_path = mat_path * dist

sum_dist_path = (np.sum(mat_path) + dist[

engineers[u], path[-1]]) / velocity + repair_time * len(path)

if sum_dist_path < sum_dist[u]:

best_path = path.copy()

creature[u] = path.copy()

if previous_best_path == best_path:

break

全部论文请见下方“ 只会建模 QQ名片” 点击QQ名片即可