2021年暑假数学建模第二次模拟赛：无人机路径规划（现代智能优化算法）

本系列赛题、数据获取：
2021年暑假数学建模模拟赛（赛题+数据+分析）
不直接提供论文等资料，分析已经很详细了
整理不易，欢迎点赞+关注+收藏

2021年暑假数学建模第二次模拟赛：无人机路径规划（现代智能优化算法）

赛题

题目细节真的很多，也很难，仔细看看吧，建议阅读详细赛题，是对阅读分析能力很好的锻炼

思路

问题一

最终结果是一个661*621的表格的大网内的通信覆盖情况，直接对每一个点进行求解。让其与这16基站进行通讯。
3种通讯方式和三种高度实际上是9个子问题，对每个子问题
1.遍历每一个空中的点
2.尝试连接基站x
3.查询基站x是否支持通信方式dx，若不支持返回2，尝试下一个
4.将该基站点转化为经纬度(问题中提供的经纬度是其高程图原点)
5.2计算距离，若距离超过最远通信距离，返回2
6.将该点映射到地面格点，空中点相当于在第a1行第b1列
7.求a0，b0到a1，b1直线
8.对于直线经过的格点，获取该点高程h0，获取该线在空中应该的高度h1
9.比较h0，h1，发现h1＜h0，说明被阻，回到2尝试下一个基站

问题二

利用遗传算法求解，详细的思路可见我的博文
遗传算法求解无人机路径多目标规划问题（python实现）

问题三

其实就是在一个图里给定位置根据一定条件再画三个圆（通信覆盖网），然后就转化为问题二了
都看到这里了，点个赞再走啊

问题一

python最大的优势就是开源，可以使用python里的geopy库进行方便的经纬度转换

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from time import time
from PIL import Image

from matplotlib.font_manager import FontProperties
from pandas_profiling import ProfileReport

import geopy.distance
from geopy.distance import lonlat, distance

ARC = 399406380 /(90 * 4)

def wtd(w):
    return 10 * np.log10(w) + 30 
    
def dis(pt,pr,gt,gr,mu,f):
    lgd = (wtd(pt) - pr +gt + gr + 20 * np.log10(mu) - 32.44 - 20 * np.log10(f))/20.0
    d = np.power(10,lgd)
    return d

def calc(lon1,lat1,alt1,lon2,lat2,alt2):
    newport_ri_xy = (lon1, lat1)
    cleveland_oh_xy = (lon2, lat2)
    gdist = distance(lonlat(*newport_ri_xy), lonlat(*cleveland_oh_xy)).m
    acw_dist = np.sqrt(gdist**2 + (alt1 - alt2)**2)
    return acw_dist

def calcxy(lon1,lat1,lon2,lat2):
    x = (lon2 - lon1) * ARC * np.cos(np.radians(lat1))
    y = (lat2 - lat1) * ARC
    return x,y

dem = [np.nan for i in range(16)]
dem[0] = pd.read_excel('dem.xlsx', header = None, sheet_name = 'A01')
dem[1] = pd.read_excel('dem.xlsx', header = None, sheet_name = 'A02')
dem[2] = pd.read_excel('dem.xlsx', header = None, sheet_name = 'A03')
dem[3] = pd.read_excel('dem.xlsx', header = None, sheet_name = 'A04')
dem[4] = pd.read_excel('dem.xlsx', header = None, sheet_name = 'A05')
dem[5] = pd.read_excel('dem.xlsx', header = None, sheet_name = 'A06')
dem[6] = pd.read_excel('dem.xlsx', header = None, sheet_name = 'B01')
dem[7] = pd.read_excel('dem.xlsx', header = None, sheet_name = 'B02')
dem[8] = pd.read_excel('dem.xlsx', header = None, sheet_name = 'B03')
dem[9] = pd.read_excel('dem.xlsx', header = None, sheet_name = 'B04')
dem[10] = pd.read_excel('dem.xlsx', header = None, sheet_name = 'B05')
dem[11] = pd.read_excel('dem.xlsx', header = None, sheet_name = 'B06')
dem[12] = pd.read_excel('dem.xlsx', header = None, sheet_name = 'B07')
dem[13] = pd.read_excel('dem.xlsx', header = None, sheet_name = 'B08')
dem[14] = pd.read_excel('dem.xlsx', header = None, sheet_name = 'B09')
dem[15] = pd.read_excel('dem.xlsx', header = None, sheet_name = 'B10')

info = pd.read_excel('info.xlsx',index_col = 0)

dsd = np.zeros(3)
dsd[0] = dis(150,-100,0,0,0.9,1624)
dsd[1] = dis(80,-90,0,0,0.9,400)
dsd[2] = dis(70,-88,0,0,0.9,360)
print('D1/D2/D3通信距离(单位为km)分别是：',dsd[0],dsd[1],dsd[2])

stime = time()
cnt = 0
gda = [3000,6000,9000]
mp = [np.zeros((661,621)) for i in range(9)]
for gd in gda: 
    for tx in range(3):
        for i in range(621):
            for j in range(661):
                x = 142 + i * 0.05
                y = 11 + j * 0.05
                for jz in range(16):
                    cnt += 1
                    if info[tx][jz] == 0:
                        continue
                        
                    rx = info['x'][jz]
                    ry = info['y'][jz]                  
                    ijzlon = int(rx/info['d'][jz])
                    ijzlat = int(ry/info['d'][jz])
                    
                    disztj = np.sqrt(rx ** 2 + ry ** 2)
                    theztj = np.degrees(np.arctan(rx/ry))
                    jzlonlat = geopy.distance.distance(meters = disztj).destination((info['B'][jz], info['L'][jz]), bearing = theztj)
                    
                    jzlon = jzlonlat[1]
                    jzlat = jzlonlat[0]
                    jzh = info['h'][jz] + dem[jz][ijzlat][ijzlon]

                    if calc(x,y,gd,jzlon,jzlat,jzh) > dsd[tx] * 1000:
                        continue
                    
                    #print(gd,tx,i,j,jz,x,y,gd,jzlon,jzlat,jzh,calc(x,y,gd,jzlon,jzlat,jzh))
                    disx,disy = calcxy(x,y,jzlon,jzlat)
                    gdx,gdy = disx / info['d'][jz], disy / info['d'][jz]
                    
                    l = int(ijzlon)
                    r = min(int(gdx+0.99999),dem[jz].shape[1])
                    if int(jzlon)>int(gdx+0.99999):
                        rg = [x for x in range(l,r,-1)]
                    else:
                        rg = [x for x in range(l,r)]
                    
                    k = (gdy - ijzlat) / (gdx - ijzlon) 
                    b = ijzlat - k * ijzlon
                    flag = True
                    for x0 in rg:
                        y0 = int(k * x0 + b)
                        if y0 >= dem[jz].shape[0]:
                            #print(0,gd,tx,i,j,jz)
                            break
                            
                        if y0 < 0:
                            #print(1,gd,tx,i,j,jz)
                            continue
                        
                        if x0 < 0:
                            continue
                            
                        hkz = calc(x,y,gd,jzlon,jzlat,jzh) * (x0 - l)/(int(gdx+0.99999) - l)
                        print(k,b,x0,y0,l,r,hkz,dem[jz][y0][x0] - jzh, dem[jz][y0][x0] ,jzh,ijzlon,ijzlat)
                        if hkz < dem[jz][y0][x0] - jzh:
                            flag = False
                            break 
                                       
                    if flag:
                        sheet = tx * 3 + int((gd - 3000)/3000)
                        mp[sheet][j][i] = 1
                        
                        if cnt % 100 == 0:
                            etime = time() - stime
                            m, s = divmod(etime, 60)
                            h, m = divmod(m, 60)
                            print('当前进度：{:.2f}%'.format(cnt * 100 / (3 * 3 * 661 * 621 * 16)))
                            print('已用时间:'+'{}:{}:{}'.format(h,m,np.around(s,decimals = 2)))
                            img = Image.fromarray(mp[sheet]*255)
                            img = img.convert('L')
                            img.save(str(sheet) + '.png')
                                      
pd.DataFrame(mp[0]).to_csv('D1_3000.csv')
pd.DataFrame(mp[1]).to_csv('D1_6000.csv')
pd.DataFrame(mp[2]).to_csv('D1_9000.csv')
pd.DataFrame(mp[3]).to_csv('D2_3000.csv')
pd.DataFrame(mp[4]).to_csv('D2_6000.csv')
pd.DataFrame(mp[5]).to_csv('D2_9000.csv')
pd.DataFrame(mp[6]).to_csv('D3_3000.csv')
pd.DataFrame(mp[7]).to_csv('D3_6000.csv')
pd.DataFrame(mp[8]).to_csv('D3_9000.csv')

img = Image.fromarray(mp[0]*255)
img = img.convert('L')
img.save('D1_3000.png')
img = Image.fromarray(mp[1]*255)
img = img.convert('L')
img.save('D1_6000.png')
img = Image.fromarray(mp[2]*255)
img = img.convert('L')
img.save('D1_9000.png')
img = Image.fromarray(mp[3]*255)
img = img.convert('L')
img.save('D2_3000.png')
img = Image.fromarray(mp[4]*255)
img = img.convert('L')
img.save('D2_6000.png')
img = Image.fromarray(mp[5]*255)
img = img.convert('L')
img.save('D2_9000.png')
img = Image.fromarray(mp[6]*255)
img = img.convert('L')
img.save('D3_3000.png')
img = Image.fromarray(mp[7]*255)
img = img.convert('L')
img.save('D3_6000.png')
img = Image.fromarray(mp[8]*255)
img = img.convert('L')
img.save('D3_9000.png')

最终结果应该是这样子的

问题二三

细节和要素过多，单独写了一篇博文
遗传算法求解无人机路径多目标规划问题（python实现）