赞
踩
本系列赛题、数据获取:
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实现)
Copyright © 2003-2013 www.wpsshop.cn 版权所有,并保留所有权利。