第一问 参考司守奎老师《python数学实验与建模》（其实就是改了改数据）
1. import numpy as np 2. import pandas as pd 3. from numpy.random import randint, rand, shuffle 4. from matplotlib.pyplot import plot, show, rc 5. from numpy import loadtxt,radians,sin,cos,inf,exp 6. from numpy import array,r_,c_,arange,savetxt 7. from numpy.lib.scimath import arccos 8. data=pd.read_excel(r"D:\数学建模\赛题\20深圳杯\C题附件2.xlsx",encoding="gbk") 9. 10. x=data["传感器经度"];y=data["传感器纬度"]; 11. d1=([[120.7015202,36.37423]]); xy=c_[x,y] 12. xy=r_[xy,d1]; N=xy.shape[0] 13. t=radians(xy) #转化为弧度 14. d=array([[6370*arccos(cos(t[i,0]-t[j,0])*cos(t[i,1])*cos(t[j,1])+ 15. sin(t[i,1])*sin(t[j,1])) for i in range(N)] 16. for j in range(N)]).real 17. savetxt('各节点间的实际距离.txt',c_[xy,d]) #把数据保存到文本文件，供下面使用 18. 19. 20. a=np.loadtxt("各节点间的实际距离.txt") 21. xy,d=a[:,:2],a[:,2:]; N=len(xy) 22. w=50; g=10 #w为种群的个数，g为进化的代数 23. J=[]; 24. 25. for i in np.arange(w): 26. c=np.arange(1,N-1); shuffle(c) 27. c1=np.r_[1,c,N-1]; flag=1 28. while flag>0: 29. flag=0 30. for m in np.arange(1,N-3): 31. for n in np.arange(m+1,N-2): 32. if d[c1[m],c1[n]]+d[c1[m+1],c1[n+1]]<d[c1[m],c1[m+1]]+d[c1[n],c1[n+1]]: 33. c1[m+1:n+1]=c1[n:m:-1]; flag=1 34. c1[c1]=np.arange(N); J.append(c1) 35. J=np.array(J)/(N-1) 36. for k in np.arange(g): 37. A=J.copy() 38. c1=np.arange(w); shuffle(c1) #交叉操作的染色体配对组 39. c2=randint(2,100,w) #交叉点的数据 40. for i in np.arange(0,w,2): 41. temp=A[c1[i],c2[i]:N-1] #保存中间变量 42. A[c1[i],c2[i]:N-1]=A[c1[i+1],c2[i]:N-1] 43. A[c1[i+1],c2[i]:N-1]=temp 44. B=A.copy() 45. by=[] #初始化变异染色体的序号 46. while len(by)<1: by=np.where(rand(w)<0.1) 47. by=by[0]; B=B[by,:] 48. G=np.r_[J,A,B] 49. ind=np.argsort(G,axis=1) #把染色体翻译成0,1，…，101 50. NN=G.shape[0]; L=np.zeros(NN) 51. for j in np.arange(NN): 52. for i in np.arange(30): 53. L[j]=L[j]+d[ind[j,i],ind[j,i+1]] 54. ind2=np.argsort(L) 55. J=G[ind2,:] 56. path=ind[ind2[0],:]; zL=L[ind2[0]] 57. xx=xy[path,0]; yy=xy[path,1]; rc('font',size=16) 58. xx_yy=c_[xx,yy] 59. plot(xx,yy,'-*'); show() #画路径 60. print("所求的路径长度为：",zL) 61. 62. 63. from scipy.integrate import quad 64. data1=pd.DataFrame(xx_yy) 65. data1.columns = ['传感器经度','传感器纬度'] 66. a=[] 67. for i in range(0,31): a.append(i) 68. data1['顺序']=a 69. 70. data2 = pd.merge(data,data1,how = 'left',on=['传感器经度','传感器纬度']) 71. data2=data2.sort_values(by='顺序') 72. row={'节点':'数据中心','传感器经度':120.701520 ,'传感器纬度':36.374227,'能量消耗速率(mA/h)':'/','顺序':30} 73. data2=data2.append(row,ignore_index=True) 74. data2.to_excel('path.xls') #保存到excel第二问
1. import numpy as np 2. import pandas as pd 3. from numpy.random import randint, rand, shuffle 4. from matplotlib.pyplot import plot, show, rc 5. from numpy import loadtxt,radians,sin,cos,inf,exp 6. from numpy import array,r_,c_,arange,savetxt 7. from numpy.lib.scimath import arccos 8. 9. data=pd.read_excel(r"D:\数学建模\赛题\20深圳杯\path.xls",encoding="gbk") 10. 11. 12. data.iloc[0,4] = 0 13. data.iloc[30,4] = 0 14. 15. seit = data['能量消耗速率(mA/h)']/3600 16. 17. from sympy import * 18. E=symbols('E') 19. # r=symbols('r') 20. v=4.167 21. r=0.5 22. f=40 23. data['ta'] = data['距离']/v #每一段路上的时间 24. 25. tb = [0 for index in range(31)] 26. tb[1] = seit[1]*data['ta'][1]/r 27. 28. def tb_(x): 29. t_1=data['ta'][1]; t_2=0 30. for i in range(2,x): 31. for n in range(1,i+1): 32. t_1 += data['ta'][n] 33. for m in range(1,i): 34. t_2 += tb[i-1] 35. tb[i]=together((seit[i]*(t_1+t_2))/r) 36. tb_(30) 37. 38. data['tb'] = tb 39. 40. # 第一次充电 41. a = [] 42. for i in range(1,30): 43. a_ = v*(sum(data['ta'][1:i])+sum(tb[0:i-1]))+f 44. a.append(a_) 45. 46. data2 = { 47. "节点":data['节点'][1:30], 48. "需要的最小电容":a} 49. data2=pd.DataFrame(data2) 50. 51. data2.to_excel('第一次充电各传感器需要的最小电容.xls') 52. 53. # 第二次充电 54. b = [] 55. for i in range(1,31): 56. b_ = v*(sum(data['ta'][1:i])+sum(tb[0:i-1]))+f 57. b.append(b_) 58. b[29] 59. 60. max(a)第三问 就费劲了 前两问做了一天半，第三问做了一周，他们俩到处给我找文献分析，撸了一套遗传算法出来
1. #导包 2. import numpy as np 3. import pandas as pd 4. from matplotlib import pyplot as plt 5. import math 6. import copy 7. import random 8. 9. from numpy.random import randint, rand, shuffle 10. from matplotlib.pyplot import plot, show, rc 11. from numpy import loadtxt,radians,sin,cos,inf,exp 12. from numpy import array,r_,c_,arange,savetxt 13. from numpy.lib.scimath import arccos 14. from math import radians, cos, sin, asin, sqrt 15. 16. data=pd.read_excel(r"C题附件2.xlsx") 17. 18. # 初始化 19. def xy(): 20. li = [] 21. data=pd.read_excel(r"C题附件2.xlsx") 22. data=pd.DataFrame(data) 23. x=data['传感器经度'];y=data['传感器纬度'] 24. for i in range(30): 25. li.append(np.array([x[i],y[i]])) 26. li.append(np.array([x[0],y[0]]));li.append(np.array([x[0],y[0]]));li.append(np.array([x[0],y[0]])) #增设三个虚点 27. li = np.array(li) 28. return li 29. li = xy() 30. 31. # 计算距离 32. def haversine(lon1, lat1, lon2, lat2): # 经度1，纬度1，经度2，纬度2 （十进制度数） 33. 34. # 将十进制度数转化为弧度 35. lon1, lat1, lon2, lat2 = map(radians, [lon1, lat1, lon2, lat2]) 36. 37. # haversine公式 38. a = sin((lat2-lat1)/2)**2 + cos(lat1) * cos(lat2) * sin((lon2-lon1)/2)**2 39. c = 2 * asin(sqrt(a)) 40. r = 6370 # 地球平均半径，单位为公里 41. return c * r#输出单位km 42. 43. def DMAT(): 44. data=pd.read_excel(r"C题附件2.xlsx") 45. D=np.zeros((len(li),len(li))) 46. for i in range(0,len(li)): 47. for j in range(0,len(li)): 48. D[i][j]=haversine(li[:,0][i],li[:,1][i],li[:,0][j],li[:,1][j]) 49. # print(D[i][j]) 50. D[0][30]=float('inf');D[0][31]=float('inf');D[0][32]=float('inf') 51. D[30][0]=float('inf');D[31][0]=float('inf');D[32][0]=float('inf') 52. D[30][31]=float('inf');D[30][32]=float('inf') 53. D[31][30]=float('inf');D[31][32]=float('inf') 54. D[32][30]=float('inf');D[32][31]=float('inf') #各虚点与数据中心之间，虚点与虚点之间距离设为无穷大 55. return D 56. # DMAT=pd.DataFrame(DMAT()) 57. DMAT = DMAT() 58. 59. li = [i for i in range(1,33)] 60. print(li[1:33]) 61. 62. # 初始化，随机排列生成100条路径，以3个0为切断点将每条路径分为4条路径,给每条路径加上起点终点 63. def population(): 64. num=100 65. population=[] 66. for i in range(num): 67. random.shuffle(li) 68. population.append(copy.deepcopy(li)) 69. for i in range(num): 70. population[i].insert(0,0) 71. population[i].append(0) 72. return population 73. population=population() 74. 75. # 适应度计算 76. def aimFunction(routes, DMAT): 77. distance = 0 78. for i in range(len(routes)-1): 79. distance += DMAT[routes[i]][routes[i+1]] 80. return 1.0/distance 81. 82. def fitness(population): 83. value=[] 84. for i in range(len(population)): 85. value.append(aimFunction(population[i],DMAT)) 86. if value[i]<0: 87. value[i]=0 88. return value 89. value = fitness(population) 90. 91. # 选择 92. def selection(population, value): # 轮盘赌选择 93. fitness_sum = [] 94. for i in range(len(value)): 95. if i == 0: 96. fitness_sum.append(value[i]) 97. else: 98. fitness_sum.append(fitness_sum[i - 1] + value[i]) 99. for i in range(len(fitness_sum)): 100. fitness_sum[i] /= sum(value) 101. # select new population 102. population_new = [] 103. for i in range(len(value)): 104. rand = np.random.uniform(0, 1) 105. for j in range(len(value)): 106. if j == 0: 107. if 0 < rand and rand <= fitness_sum[j]: 108. population_new.append(population[j]) 109. else: 110. if fitness_sum[j - 1] < rand and rand <= fitness_sum[j]: 111. population_new.append(population[j]) 112. return population_new 113. population_new = selection(population, value) 114. 115. 116. # 交叉 117. def amend(entity, low, high): 118. """ 119. 修正个体 120. :param entity: 个体 121. :param low: 交叉点最低处 122. :param high: 交叉点最高处 123. :return: 124. """ 125. length = len(entity) 126. cross_gene = entity[low:high] # 交叉基因 127. not_in_cross = [] # 应交叉基因 128. raw = entity[0:low] + entity[high:] # 非交叉基因 129. 130. for i in range(1,length+1): 131. if not i in cross_gene: 132. not_in_cross.append(i) 133. 134. error_index = [] 135. for i in range(len(raw)): 136. if raw[i] in not_in_cross: 137. not_in_cross.remove(raw[i]) 138. else: 139. error_index.append(i) 140. for i in range(len(error_index)): 141. raw[error_index[i]] = not_in_cross[i] 142. 143. entity = raw[0:low] + cross_gene + raw[low:] 144. return entity 145. 146. 147. 148. def crossover(population_new, pc): 149. """ 150. 交叉 151. :param population_new: 种群 152. :param pc: 交叉概率 153. :return: 154. """ 155. half = int(len(population_new) / 2) 156. father = population_new[:half] 157. mother = population_new[half:] 158. np.random.shuffle(father) 159. np.random.shuffle(mother) 160. offspring = [] 161. for i in range(half): 162. if np.random.uniform(0, 1) <= pc: 163. cut1 = 0 164. cut2 = np.random.randint(0, len(population_new[0])) 165. if cut1 > cut2: 166. cut1, cut2 = cut2, cut1 167. if cut1 == cut2: 168. son = father[i] 169. daughter = mother[i] 170. else: 171. son = father[i][0:cut1] + mother[i][cut1:cut2] + father[i][cut2:] 172. son = [0]+amend(son[1:(len(son)-1)], cut1, cut2)+[0] 173. daughter = mother[i][0:cut1] + father[i][cut1:cut2] + mother[i][cut2:] 174. daughter = [0]+amend(daughter[1:(len(daughter)-1)], cut1, cut2)+[0] 175. else: 176. son = father[i] 177. daughter = mother[i] 178. offspring.append(son) 179. offspring.append(daughter) 180. return offspring 181. 182. 183. # 变异 184. def mutation(offspring, pm): 185. for i in range(len(offspring)): 186. if np.random.uniform(0, 1) <= pm: 187. position1 = np.random.randint(1, len(offspring[i])-1) 188. position2 = np.random.randint(1, len(offspring[i])-1) 189. offspring[i][position1],offspring[i][position2] = offspring[i][position2],offspring[i][position1] 190. return offspring 191. 192. 193. import numpy as np 194. 195. 196. #主逻辑代码 197. def show_population(population): 198. # x = [i / 100 for i in range(900)] 199. x = [i / 100 for i in range(-450, 450)] 200. y = [0 for i in range(900)] 201. for i in range(900): 202. y[i] = aimFunction(x[i]) 203. 204. population_10 = [decode(p) for p in population] 205. y_population = [aimFunction(p) for p in population_10] 206. 207. plt.plot(x, y) 208. plt.plot(population_10, y_population, 'ro') 209. plt.show() 210. 211. 212. if __name__ == '__main__': 213. 214. locations = xy() 215. 216. t = [] 217. route = [] 218. for i in range(4001): 219. # selection 220. value = fitness(population) 221. population_new = selection(population, value) 222. # crossover 223. offspring = crossover(population_new, 0.65) 224. # mutation 225. population = mutation(offspring, 0.02) 226. result = [] 227. for j in range(len(population)): 228. result.append(1.0 / aimFunction(population[j], DMAT)) 229. route.append(population[j]) 230. t.append(min(result)) 231. if i % 400 == 0: 232. min_entity = population[result.index(min(result))] 233. print(min_entity) 234. plt.show() 235. plt.plot(t) 236. plt.show() 237. 238. 239. # 总距离 240. route_f=[0, 20, 18, 7, 9, 8, 15, 11, 14, 6, 16, 27, 13, 3, 31, 4, 32, 5, 10, 12, 2, 1, 30, 17, 29, 26, 25, 19, 21, 23, 28, 24, 22, 0] 241. distance = 1.0/aimFunction(route_f, DMAT) 242. distance 243. 244. 245. distance_e = [] 246. for i in range(1,len(route_f)): 247. if i>0: 248. distance_e.append(DMAT[route_f[i-1]][route_f[i]]) 249. distance_e.insert(0,0) 250. 251. n = [] 252. for i in range(30): 253. n.append(i) 254. data['节点']=n 255. 256. x = [] 257. for i in range(34): 258. x.append(i) 259. 260. # data1 = pd.DataFrame({'节点':pd.Series(route_f),'距离':pd.Series(distance_e),'顺序':pd.Series(x)}) 261. data1 = pd.DataFrame({'节点':(route_f),'距离':(distance_e),'顺序':(x)}) 262. 263. data1 = pd.merge(data1,data,how = 'left',on='节点') 264. data1 = data1.sort_values(by='顺序',ascending=True) 265. 266. data1['传感器经度'].fillna(120.7015202, inplace=True) 267. data1['传感器纬度'].fillna(36.37423 , inplace=True) 268. 269. 270. from sympy import * 271. 272. E=symbols('E') 273. # r=symbols('r') 274. v=4.167 275. r=0.5 276. f=40 277. 278. data1['ta'] = data1['距离']/v #每一段路上的时间 279. 280. data1['能量消耗速率(mA/h)'].fillna(0, inplace=True) 281. seit = list(data1['能量消耗速率(mA/h)']) 282. seit[0]=0; 283. seit[33]=0 284. 285. for i in range(len(seit)): 286. seit[i]/=3600 287. 288. data_1 = data1.loc[0:14] 289. data_1['节点'][14]=0 290. 291. # 充电时间1 292. tb = [0 for index in range(len(data_1))] 293. tb[1] = seit[1]*data_1['ta'][1]/r 294. tb[1] 295. 296. def tb_(x): 297. t_1=data_1['ta'][1]; t_2=0 298. for i in range(2,x): 299. for n in range(1,i+1): 300. t_1 += data_1['ta'][n] 301. for m in range(1,i): 302. t_2 += tb[i-1] 303. tb[i]=together((seit[i]*(t_1+t_2))/r) 304. tb_(len(data_1)) 305. data_1['tb'] = tb 306. data_1['节点'][14]=0 307. data_1.to_excel('第一辆移动充电器路径资料.xls') 308. 309. # 第一次充电 第一辆车 310. a = [] 311. for i in range(1,len(data_1)-1): 312. a_ = E-v*(sum(data_1['ta'][1:i])+sum(tb[0:i-1]))-f 313. a.append(a_) 314. 315. # 第二次充电 第一辆车 316. b = [] 317. for i in range(1,len(data_1)): 318. b_ = E-v*(sum(data_1['ta'][1:i])+sum(tb[0:i-1]))-f 319. b.append(b_) 320. b[len(data_1)-2] 321. 322. data_2=data1.loc[14:16] 323. 324. # 充电时间2 325. tb = [0 for index in range(len(data_2))] 326. tb[1] = seit[15]*0.099537/r 327. tb[1] 328. data_2['tb'] = tb 329. data_2['节点'][0]=0 330. data_2['节点'][2]=0 331. data_2.to_excel('第二辆移动充电器路径资料.xls') 332. 333. # 第一次充电 第二辆车 334. a = [] 335. for i in range(1,len(data_2)-1): 336. a_ = E-v*(sum(data_2['ta'][1:i])+sum(tb[0:i-1]))-f 337. a.append(a_) 338. 339. # 第二次充电 第二辆车 340. b = [] 341. for i in range(1,len(data_2)): 342. b_ = E-v*(sum(data_2['ta'][1:i])+sum(tb[0:i-1]))-f 343. b.append(b_) 344. b[len(data_2)-2] 345. 346. data_3 = data1.loc[16:22] 347. 348. # 充电时间3 349. seit3=seit[16:23] 350. tb = [0 for index in range(len(data_3))] 351. tb[1] = seit3[1]*list(data_3['ta'])[1]/r 352. tb[1] 353. 354. def tb_(x): 355. t_1=list(data_3['ta'])[1]; t_2=0 356. for i in range(2,x): 357. for n in range(1,i+1): 358. t_1 += list(data_3['ta'])[n] 359. for m in range(1,i): 360. t_2 += tb[i-1] 361. tb[i]=together((seit3[i]*(t_1+t_2))/r) 362. tb_(len(data_3)) 363. data_3['tb'] = tb 364. data_3['节点'][16]=0 365. data_3['节点'][22]=0 366. data_3.to_excel('第三辆移动充电器路径资料.xls') 367. 368. # 第一次充电 第三辆车 369. a = [] 370. for i in range(1,len(data_3)-1): 371. a_ = E-v*(sum(data_3['ta'][1:i])+sum(tb[0:i-1]))-f 372. a.append(a_) 373. 374. # 第二次充电 第三辆车 375. b = [] 376. for i in range(1,len(data_3)): 377. b_ = E-v*(sum(data_3['ta'][1:i])+sum(tb[0:i-1]))-f 378. b.append(b_) 379. b[len(data_3)-2] 380. 381. data_4=data1.loc[22:33] 382. 383. # 充电时间4 384. seit4=seit[22:] 385. tb = [0 for index in range(len(data_4))] 386. tb[1] = seit4[1]*list(data_4['ta'])[1]/r 387. tb[1] 388. 389. def tb_(x): 390. t_1=list(data_4['ta'])[1]; t_2=0 391. for i in range(2,x): 392. for n in range(1,i+1): 393. t_1 += list(data_4['ta'])[n] 394. for m in range(1,i): 395. t_2 += tb[i-1] 396. tb[i]=together((seit4[i]*(t_1+t_2))/r) 397. tb_(len(data_4)) 398. data_4['tb'] = tb 399. data_4['节点'][22]=0 400. data_4.to_excel('第四辆移动充电器路径资料.xls') 401. 402. # 第一次充电 第四辆车 403. a = [] 404. for i in range(1,len(data_4)-1): 405. a_ = E-v*(sum(data_4['ta'][1:i])+sum(tb[0:i-1]))-f 406. a.append(a_) 407. 408. # 第二次充电 第四辆车 409. b = [] 410. for i in range(1,len(data_4)): 411. b_ = E-v*(sum(data_4['ta'][1:i])+sum(tb[0:i-1]))-f 412. b.append(b_) 413. b[len(data_4)-2] 414. 415. #路径图 416. x1=data_1['传感器经度'];y1=data_1['传感器纬度'] 417. x2=data_2['传感器经度'];y2=data_2['传感器纬度'] 418. x3=data_3['传感器经度'];y3=data_3['传感器纬度'] 419. x4=data_4['传感器经度'];y4=data_4['传感器纬度'] 420. 421. plot(x1,y1,'-*',color='lightcoral'); plot(x2,y2,'-*',color='green'); 422. plot(x3,y3,'-*',color='gold'); plot(x4,y4,'-*',color='royalblue'); 423. show() #画路径有时间再来写思路
