UAV-RIS-MEC附python代码.zip

import math import numpy as np import math as mt M = 100 phase_shift = 100 B = 8000 # 8MHz Xi = mt.pow(10, (3 / 10)) white_noise = 10 Power_ue = 3 * mt.pow(10, 3) # 300mW Power_uav = 5 * mt.pow(10, 3) # 500mW N_0 = mt.pow(10, (-169 / 10)) * (0.1 * mt.pow(10, 3)) / B omega1 = 0.05 omega2 = 1 chi = 1e-26 # 起飞 def launching(self): uav_index = -1 for uav in self.uav_cluster: uav_index += 1 # 让无人机高度抬高90m for _ in np.arange(1, 21, 1): theta = 0 # 单位时间内的飞行仰角（ 0 表示垂直起飞） fai = 0 # 单位时间内飞行轨迹与x轴正方向的水平夹角，垂直起降时 fai 取值随意 velocity = 3 # 飞行速度，m/s uav.energy = uav.energy - get_power_cost_per_meter(abs(velocity)) * abs(velocity) * 1 uav.flight_trace(theta, fai, velocity, 1) # 直线飞行 self.uav_cluster_x[uav_index].append(uav.x) self.uav_cluster_y[uav_index].append(uav.y) self.uav_cluster_z[uav_index].append(uav.z) return self.uav_cluster # 降落 def landing(self, uav, uav_index): for _ in np.arange(1, 3000, 1): theta = math.pi # 单位时间内的飞行仰角（ pi 表示垂直降落） fai = 0 # 单位时间内飞行轨迹与x轴正方向的水平夹角，垂直起降时 fai 取值随意 velocity = 3 # 飞行速度，m/s if uav.z >= 1 and (uav.z - velocity) < 1: velocity = uav.z - 1 uav.energy = uav.energy - get_power_cost_per_meter(abs(velocity)) * abs(velocity) * 1 uav.flight_trace(theta, fai, velocity, 1) # 直线飞行 self.uav_cluster_x[uav_index].append(uav.x) self.uav_cluster_y[uav_index].append(uav.y) self.uav_cluster_z[uav_index].append(uav.z) if uav.z == 1: break # 得到多点连线形成的多边形的中心点 # 按该方法得到的中心点物理意义不明，已弃用 # reference - https://www.csdn.net/tags/NtzaQg2sNjUxOTUtYmxvZwO0O0OO0O0O.html def get_center_location(clusters_points): x = 0 # lon y = 0 # lat z = 0 length = len(clusters_points) for lon, lat in clusters_points: lon = math.radians(float(lon)) # radians(float(lon)) Convert angle x from degrees to radians # 把角度 x 从度数转化为 弧度 lat = math.radians(float(lat)) x += math.cos(lat) * math.cos(lon) y += math.cos(lat) * math.sin(lon) z += math.sin(lat) x = float(x / length) y = float(y / length) z = float(z / length) return math.degrees(math.atan2(y, x)), math.degrees(math.atan2(z, math.sqrt(x * x + y * y))) # Gauss's area formula 高斯面积计算，多边形形心计算中使用 def cal_area(vertices): # Gauss's area formula 高斯面积计算 A = 0.0 point_p = vertices[-1] for point in vertices: A += (point[1]*point_p[0] - point[0]*point_p[1]) point_p = point return abs(A)/2 # 得到多点连线形成的多边形的形心，需要逆时针输入各点坐标，输入格式：[[-1., -1.], [-2., -1.], [-2., -2.], [-1., -2.]] # reference - https://blog.csdn.net/kindlekin/article/details/121318530 def cal_centroid(points): A = cal_area(points) c_x, c_y = 0.0, 0.0 point_p = points[-1] # point_p 表示前一节点 for point in points: c_x += ((point[0] + point_p[0]) * (point[1]*point_p[0] - point_p[1]*point[0])) c_y += ((point[1] + point_p[1]) * (point[1]*point_p[0] - point_p[1]*point[0])) point_p = point return c_x / (6*A), c_y / (6*A) # 根据已知两点坐标，求过这两点的直线解析方程： a1*x+b1*y+c1 = 0 (a >= 0) & a2*x+b2*z+c2 = 0 (a >= 0) def get_linear_equation(ue, uav): [p1x, p1y, p1z, p2x, p2y, p2z] = [ue.x, ue.y, ue.z, uav.x, uav.y, uav.z] sign1 = 1 a1 = p2y - p1y if a1 < 0: sign1 = -1 a1 = sign1 * a1 b1 = sign1 * (p1x - p2x) c1 = sign1 * (p1y * p2x - p1x * p2y) sign2 = 1 a2 = p2z - p1z if a2 < 0: sign2 = -1 a2 = sign2 * a2 b2 = sign2 * (p1x - p2x) c2 = sign2 * (p1z * p2x - p1x * p2z) return [a1, b1, c1, a2, b2, c2] # 视距判断 def line_of_sight_judgement(ue_cluster, uav_cluster, building_cluster): los_judgement = np.ones((len(ue_cluster), len(uav_cluster))) rows_index = -1 for rows in los_judgement: cols_index = -1 rows_index += 1 for cols in rows: cols_index += 1 coefficients = get_linear_equation(ue_cluster[rows_index], uav_cluster[cols_index]) for x_sample in np.arange(ue_cluster[rows_index].x, uav_cluster[cols_index].x, 0.1): y_sample = (- coefficients[2] - coefficients[0] * x_sample) / coefficients[1] z_sample = (- coefficients[5] - coefficients[3] * x_sample) / coefficients[4] for building in building_cluster: if building.x <= x_sample <= (building.x + building.dx) and \ building.y <= y_sample <= (building.y + building.dy) and \ building.z <= z_sample <= (building.z + building.dz): los_judgement[rows_index][cols_index] = 0 return los_judgement # 坠机判断（无人机与建筑物相撞或无人机之间相撞） def uav_crash_judgement(uav_cluster, building_cluster): crash_flag1 = False uav_uav_distance = get_clusters_distance(uav_cluster, uav_cluster) n = len(uav_cluster) # 无人机不会与其自身碰撞，直接给一个较大的距离 for i in range(n): uav_uav_distance[i][i] = 1000 for uav in uav_cluster: for building in building_cluster: if building.x <= uav.x <= (building.x + building.dx) and \ building.y <= uav.y <= (building.y + building.dy) and \ building.z <= uav.z <= (building.z + building.dz): crash_flag1 = True crash_flag2 = (uav_uav_distance <= 1).any() crash_judgement = crash_flag1 or crash_flag2 return crash_judgement # agent_cluster1 与 agent_cluster2 之间的距离 def get_clusters_distance(agent_cluster1, agent_cluster2): # https://blog.csdn.net/Tan_HandSome/article/details/82501902 agent_cluster1_locations = [] agent_cluster2_locations = [] for agent1 in agent_cluster1: agent_cluster1_locations.append(agent1.xyz) for agent2 in agent_cluster2: agent_cluster2_locations.append(agent2.xyz) A = np.array(agent_cluster1_locations) B = np.array(agent_cluster2_locations) BT = B.transpose() vecProduct = np.dot(A, BT) # dot production SqA = A ** 2 # square of every element in A sumSqA = np.matrix(np.sum(SqA, axis=1)) sumSqAEx = np.tile(sumSqA.transpose(), (1, vecProduct.shape[1])) SqB = B ** 2 sumSqB = np.sum(SqB, axis=1) sumSqBEx = np.tile(sumSqB, (vecProduct.shape[0], 1)) SqED = sumSqBEx + sumSqAEx - 2 * vecProduct SqED[SqED < 0] = 0.0 ED = np.sqrt(SqED) distance = np.array(ED) return distance # ue 到 uav 的路径损耗的计算 # 根据 ADVISOR-007 论文中的公式（1）—（6）建模计算 def get_ue_to_uav_path_loss(ue_cluster, uav_cluster, building_cluster): fc = 2e9 # 单位 HZ velocity_c = 3e8 # 光速，单位 m/s ue_uav_distance = get_clusters_distance(ue_cluster, uav_cluster) path_loss = 20 * np.log10(4 * math.pi * fc * ue_uav_distance / velocity_c) los_judgement = line_of_sight_judgement(ue_cluster, uav_cluster, building_cluster) rows_index = -1 for rows in los_judgement: cols_index = -1 rows_index += 1 for cols in rows: cols_index += 1 if cols == 0: path_loss[rows_index][cols_index] = path_loss[rows_index][cols_index] + 5 elif cols == 1: path_loss[rows_index][cols_index] = p