无人机辅助移动边缘计算的计算卸载优化:一种深度确定性策略梯度方法（6）——代码实现

参考连接：

[1] Wang Y , Fang W , Ding Y , et al. Computation offloading optimization for UAV-assisted mobile edge computing: a deep deterministic policy gradient approach[J]. Wireless Networks, 2021:1-16.doi：https://doi.org/10.1007/s11276-021-02632-z

https://morvanzhou.github.io/tutorials/

https://github.com/fangvv/UAV-DDPG/blob/main/DDPG/DDPG_without_behavior_noise/ddpg_algo.py

https://blog.csdn.net/weixin_43835470/article/details/120881273

6 代码实现

本部分主要包括 Actor Critc，DDPG，DQN，Edge_only，Local_only 算法的基本实现思路，具体细节不予展示。

6.1 问题环境

问题环境主要包括了环境的各种状态 s 、无人机 uav 相关参数、用户 ue 相关参数和一些影响变量。其主要功能是模拟整个实验环境。主要动作包括了无人机飞行、无人机运算、用户移动和用户本地运算。

6.1.1 参数设计

参数包括了场地大小，无人机位置，带宽，噪声功率等参数，包括模拟实验中的所有参数。

#################### uav ####################
height = ground_length = ground_width = 100  # 场地长宽均为100m，UAV飞行高度也是
sum_task_size = 60 * 1048576  # 总计算任务60 Mbits
loc_uav = [50, 50]  #无人机的初始位置bandwidth_nums = 1
# 1MHz = 10^6Hz
B = bandwidth_nums * 10 ** 6  # 带宽1MHz
p_noisy_los = 10 ** (-13)  # 噪声功率-100dBm    ****
p_noisy_nlos = 10 ** (-11)  # 噪声功率-80dBm    ****
flight_speed = 50.  # 飞行速度50m/s
f_ue = 6e8  # UE的计算频率0.6GHz
f_uav = 12e8  # UAV的计算频率1.2GHz
r = 10 ** (-27)  # 芯片结构对cpu处理的影响因子
s = 1000  # 单位bit处理所需cpu圈数1000
p_uplink = 0.1  # 上行链路传输功率0.1W
# alpha0 = -30  # 距离为1m时的参考信道增益-30dB
alpha0 = 1e-5  # 距离为1m时的参考信道增益-50dB = 1e-5   ****
T = 200  # 周期200s
delta_t = 5  # 1s飞行, 后4s用于悬停计算
slot_num = int(T / delta_t)  # 40个间隔
m_uav = 9.65  # uav质量/kg
e_battery_uav = 500000  # uav电池电量: 500kJ. ref: Mobile Edge Computing via a UAV-Mounted Cloudlet: Optimization of Bit Allocation and Path Planning#################### ues ####################
M = 4  # UE数量
block_flag_list = np.random.randint(0, 2, M)  # 4个ue，ue的遮挡情况
loc_ue_list = np.random.randint(0, 101, size=[M, 2])  # 位置信息:x在0-100随机
# task_list = np.random.randint(1048576, 2097153, M)    # 随机计算任务1~2Mbits
task_list = np.random.randint(1572864, 2097153, M)  # 随机计算任务1.5~2Mbits
# ue位置转移概率
# 0:位置不变; 1:x+1,y; 2:x,y+1; 3:x-1,y; 4:x,y-1
# 60%概率原地不动，10%概率向上下左右
loc_ue_trans_pro = np.array([[.6, .1, .1, .1, .1],[.6, .1, .1, .1, .1],[.6, .1, .1, .1, .1],[.6, .1, .1, .1, .1]])action_bound = [-1, 1]  # 对应tahn激活函数
action_dim = 4  # 第一位表示服务的ue id;中间两位表示飞行角度和距离；后1位表示目前服务于UE的卸载率
state_dim = 4 + M * 4  # uav 剩余电量, uav 位置, 剩余总任务大小, 所有ue 位置, 所有ue 任务大小, 所有ue 遮挡情况

6.1.2 初始化部分

初始化包括 uav 剩余电量, uav 位置, 剩余总任务大小, 所有 ue 位置, 所有 ue 任务大小, 所有 ue 遮挡情况。

def __init__(self):'''初始化环境包括uav 剩余电量, uav 位置, 剩余总任务大小, 所有ue 位置, 所有ue 任务大小, 所有ue 遮挡情况Returns-------None.'''# uav battery remain, uav loc, remaining sum task size, all ue loc, all ue task size, all ue block_flagself.start_state = np.append(self.e_battery_uav, self.loc_uav)self.start_state = np.append(self.start_state, self.sum_task_size)self.start_state = np.append(self.start_state, np.ravel(self.loc_ue_list))self.start_state = np.append(self.start_state, self.task_list)self.start_state = np.append(self.start_state, self.block_flag_list)self.state = self.start_state

6.1.3 重置

由于实验是划分时间段进行的，涉及到每个步骤执行后的重置。

def reset(self):'''重置环境Returns-------所有的环境信息'''self.reset_env()# uav battery remain, uav loc, remaining sum task size, all ue loc, all ue task size, all ue block_flagself.state = np.append(self.e_battery_uav, self.loc_uav)self.state = np.append(self.state, self.sum_task_size)self.state = np.append(self.state, np.ravel(self.loc_ue_list))self.state = np.append(self.state, self.task_list)self.state = np.append(self.state, self.block_flag_list)return self._get_obs()def reset_env(self):'''将所有的值设置为默认值Returns-------None.'''self.sum_task_size = 100 * 1048576  # 总计算任务60 Mbitsself.e_battery_uav = 500000  # uav电池电量: 500kJself.loc_uav = [50, 50]self.loc_ue_list = np.random.randint(0, 101, size=[self.M, 2])  # 位置信息:x在0-100随机self.reset_step()def reset_step(self):'''随机初始化ue 的计算任务和遮挡情况Returns-------None.'''self.task_list = np.random.randint(2621440, 3145729, self.M)  # 随机计算任务1.5~2Mbits -> 1.5~2 2~2.5 2.5~3 3~3.5 3.5~4self.block_flag_list = np.random.randint(0, 2, self.M)  # 4个ue，ue的遮挡情况

6.1.4 每段时间动作

在每个时间片段中进行的动作包括了选择用户，无人机飞行，计算任务，算出最大处理时延几个部分。

选择用户的时候可以随机选择用户，也可以类似DQN中根据对应的动作选择用户。

无人机飞行时需要计算飞行后的位置以及飞行时的能量消耗。

计算任务的时候需要考虑特殊情况，如计算任务已经完成，那么实验立刻终止，不进行后续步骤；最后一步的计算任务总量与我们规定的总任务量不一致，那么重置任务大小使其满足总计算任务量；无人机飞行时飞出场地，需要重新进行该步然后进行惩罚；以及无人机的电量不足以飞行或者计算，那么就将无人机所有电量都用完后，再进行下一片段。

计算最大处理时延后需要更新我们的奖惩，同时考虑用户移动。

    def step(self):  # 0: 选择服务的ue编号 ; 1: 方向theta; 2: 距离d; 3: offloading ratiostep_redo = False  #是否重做该步is_terminal = False  #是否终止ue_id = np.random.randint(0, self.M) #随机选择一个用户theta = 0  # 角度offloading_ratio = 0  # ue卸载率task_size = self.task_list[ue_id]  #获取该用户的计算任务block_flag = self.block_flag_list[ue_id]  #获取该用户的遮挡情况# 飞行距离dis_fly = 0  # 1s飞行距离# 飞行能耗# delta_t：5s（1s飞行，4s计算）  m_uav：uav质量e_fly = (dis_fly / (self.delta_t * 0.5)) ** 2 * self.m_uav * (self.delta_t * 0.5) * 0.5  # ref: Mobile Edge Computing via a UAV-Mounted Cloudlet: Optimization of Bit Allocation and Path Planning# ref：通过无人机装载的云计算的移动边缘计算:位分配和路径规划的优化# UAV飞行后的位置dx_uav = dis_fly * math.cos(theta)dy_uav = dis_fly * math.sin(theta)loc_uav_after_fly_x = self.loc_uav[0] + dx_uavloc_uav_after_fly_y = self.loc_uav[1] + dy_uav# 服务器计算耗能# f_uav:UAV的计算频率  s：单位bit处理所需cpu圈数t_server = offloading_ratio * task_size / (self.f_uav / self.s)  # 在UAV边缘服务器上计算时延# r:芯片结构对cpu处理的影响因子e_server = self.r * self.f_uav ** 3 * t_server  # 在UAV边缘服务器上计算耗能# 计算任务全部完成if self.sum_task_size == 0:  is_terminal = True# file_name = 'output.txt'# with open(file_name, 'a') as file_obj:#     file_obj.write("\n======== This episode is done ========")  # 本episode结束reward = 0# 最后一步计算任务和ue的计算任务不匹配elif self.sum_task_size - self.task_list[ue_id] < 0:  # 将最后一步计算任务改成总计算任务的剩余量self.task_list = np.ones(self.M) * self.sum_task_sizereward = 0step_redo = True# uav位置不对elif loc_uav_after_fly_x < 0 or loc_uav_after_fly_x > self.ground_width \or loc_uav_after_fly_y < 0 or loc_uav_after_fly_y > self.ground_length:reward = -100step_redo = True# uav电量不能支持飞行elif self.e_battery_uav < e_fly:  reward = -100# uav电量不能支持计算elif self.e_battery_uav - e_fly < e_server: reward = -100# 电量支持飞行,且计算任务合理,且计算任务能在剩余电量内计算else:  delay = self.com_delay(self.loc_ue_list[ue_id], np.array([loc_uav_after_fly_x, loc_uav_after_fly_y]),offloading_ratio, task_size, block_flag)  # 计算delayreward = delay  #reward设置为时延# 更新下一时刻状态self.e_battery_uav = self.e_battery_uav - e_fly - e_server  # uav 剩余电量self.sum_task_size -= self.task_list[ue_id]  # 剩余任务量for i in range(self.M):  # ue随机移动tmp = np.random.rand()if 0.6 < tmp <= 0.7:self.loc_ue_list[i] += [0, 1]elif 0.7 < tmp <= 0.8:self.loc_ue_list[i] += [1, 0]elif 0.8 < tmp <= 0.9:self.loc_ue_list[i] += [0, -1]elif 0.9 < tmp <= 1:self.loc_ue_list[i] += [-1, 0]else:self.loc_ue_list[i] += [0, 0]# np.clip是一个截取函数，用于截取数组中小于或者大于某值的部分，并使得被截取部分等于固定值。# np.clip(a, a_min, a_max, out=None):# a：输入矩阵；a_min：被限定的最小值，所有比a_min小的数都会强制变为a_min；# a_max：被限定的最大值，所有比a_max大的数都会强制变为a_max；out：可以指定输出矩阵的对象，shape与a相同# 限定了用户的位置np.clip(self.loc_ue_list[i], 0, 100)# self.task_list = np.random.randint(1048576, 2097153, self.M)  # ue随机计算任务1~2Mbits# 随机初始化ue 的计算任务和遮挡情况self.reset_step()return reward, is_terminal, step_redo

6.1.5 计算最大处理时延

根据用户位置，无人机位置，ue卸载率，ue的计算任务，ue的遮挡情况这些参数来计算出最大处理时延。

参照下列公式：
max⁡{tlocal,k(i),tUAV,k(i)+ttr,k(i)}\max \left\{t_{l o c a l, k}(i), t_{U A V, k}(i)+t_{t r, k}(i)\right\} max{tlocal,k(i),tUAV,k(i)+ttr,k(i)}

    def com_delay(self, loc_ue, loc_uav, offloading_ratio, task_size, block_flag):'''计算花费Parameters----------loc_ue : 用户位置loc_uav : 无人机位置offloading_ratio : ue卸载率task_size : ue的计算任务block_flag : ue的遮挡情况Returns-------时延'''# 获取uav与ue之间的距离dx = loc_uav[0] - loc_ue[0]dy = loc_uav[1] - loc_ue[1]dh = self.heightdist_uav_ue = np.sqrt(dx * dx + dy * dy + dh * dh)# 噪声功率p_noise = self.p_noisy_losif block_flag == 1:p_noise = self.p_noisy_nlosg_uav_ue = abs(self.alpha0 / dist_uav_ue ** 2)  # 信道增益trans_rate = self.B * math.log2(1 + self.p_uplink * g_uav_ue / p_noise)  # 上行链路传输速率bpst_tr = offloading_ratio * task_size / trans_rate  # 上传时延,1B=8bitt_edge_com = offloading_ratio * task_size / (self.f_uav / self.s)  # 在UAV边缘服务器上计算时延t_local_com = (1 - offloading_ratio) * task_size / (self.f_ue / self.s)  # 本地计算时延# 比较上传到服务器的时间和本地计算的时间，时延越长，返回的值越高return max([t_tr + t_edge_com, t_local_com])

6.2 状态标准化

使用状态归一化算法对观测状态进行预处理，从而更有效地训练 DNN 。该算法将每个变量的最大值与最小值之差作为尺度因子。所提出的状态归一化算法可以很好地解决输入变量的大小差异问题。

在我们的工作中，变量 Ebattery (i),q(i),p1(i),…,pK(i),Dremain (i),D1(i),…,DK−1(i)E_{\text {battery }}(i), \mathbf{q}(i), \mathbf{p}_{1}(i), \ldots, \mathbf{p}_{K}(i), D_{\text {remain }}(i), D_{1}(i), \ldots, D_{K-1}(i)Ebattery (i),q(i),p1(i),…,pK(i),Dremain (i),D1(i),…,DK−1(i) 和 DK−1(i)D_{K-1}(i)DK−1(i) 在状态集中处于不同的序列，这可能导致在训练中出现问题。如算法 1 所示，通过状态归一化对这些变量进行归一化，以防止出现这种问题。在状态归一化算法中，我们使用了五个尺度因子。每个因素可以解释如下。利用缩放因子 γb\gamma_bγb 来缩小无人机电池容量。由于 UAV 和 UE 具有相同的 x 和 y 坐标范围，我们使用 γx\gamma_xγx 和 γy\gamma_yγy 分别缩小UAV和UE的x和y坐标。我们使用 γDrm\gamma_{D_{rm}}γDrm 来缩小整个时间段内剩余的任务，使用 γDUE\gamma_{D_{UE}}γDUE 来缩小时间段 i 内每个终端的任务大小。

import numpy as np
from UAV_env import UAVEnvenv = UAVEnv()
M = env.Mclass StateNormalization(object):'''状态标准化类'''def __init__(self):self.high_state = np.array([5e5, env.ground_length, env.ground_width, 100 * 1048576])self.high_state = np.append(self.high_state, np.ones(M * 2) * env.ground_length)self.high_state = np.append(self.high_state, np.ones(M) * 3145728)self.high_state = np.append(self.high_state, np.ones(M))self.low_state = np.zeros(20)  # uav loc, ue loc, task size, block_flagself.low_state[len(self.low_state) - 2 * M:len(self.low_state) - M] = np.ones(M) * 2621440def state_normal(self, state):state[len(state) - 2 * M: len(state) - M] -= 2621440res = state / (self.high_state - self.low_state)return res

6.3 强化学习类

这里使用自己想要的强化学习类来完成对应的策略选择。以DQN为例：

# Deep Q Network off-policy
class DeepQNetwork:def __init__(self,n_actions,n_features,learning_rate=0.1,reward_decay=0.001,e_greedy=0.99,replace_target_iter=200,memory_size=MEMORY_CAPACITY,batch_size=BATCH_SIZE,# e_greedy_increment=8.684615e-05,output_graph=False,):'''初始化DQN网络Parameters----------n_actions : TYPE行为空间大小n_features : TYPE环境特征learning_rate : TYPE, optional学习率. The default is 0.1.reward_decay : TYPE, optional奖励衰减因子. The default is 0.001.e_greedy : TYPE, optionale贪心因子. The default is 0.99.replace_target_iter : TYPE, optional每多少次更新target参数. The default is 200.memory_size : TYPE, optional记忆库大小. The default is MEMORY_CAPACITY.batch_size : TYPE, optional批处理大小. The default is BATCH_SIZE.# e_greedy_increment : TYPE, optionale贪心因子增长幅度. The default is 8.684615e-05.output_graph : TYPE, optional是否输出图表. The default is False.Returns-------None.'''self.n_actions = n_actionsself.n_features = n_featuresself.lr = learning_rateself.gamma = reward_decayself.epsilon_max = e_greedyself.replace_target_iter = replace_target_iterself.memory_size = memory_sizeself.batch_size = batch_size# self.epsilon_increment = e_greedy_increment# self.epsilon = 0 if e_greedy_increment is not None else self.epsilon_maxself.epsilon = 0.99# total learning stepself.learn_step_counter = 0# initialize zero memory [s, a, r, s_]# 初始化记忆库self.memory = np.zeros((MEMORY_CAPACITY, s_dim * 2 + 2), dtype=np.float32)  # memory里存放当前和下一个state，动作和奖励# consist of [target_net, evaluate_net]# 构建神经网络self._build_net()# tf.get_collection()：从一个集合中取出变量 # tf.GraphKeys 包含所有graph collection中的标准集合名# tf.GraphKeys.GLOBAL_VARIABLES 则应该是所有的图变量的集合t_params = tf.get_collection(tf.GraphKeys.GLOBAL_VARIABLES, scope='target_net')e_params = tf.get_collection(tf.GraphKeys.GLOBAL_VARIABLES, scope='eval_net')with tf.variable_scope('hard_replacement'):self.target_replace_op = [tf.assign(t, e) for t, e in zip(t_params, e_params)]self.sess = tf.Session()if output_graph:# $ tensorboard --logdir=logstf.summary.FileWriter("logs/", self.sess.graph)# 全局变量初始化self.sess.run(tf.global_variables_initializer())self.cost_his = []def _build_net(self):'''构建所有的网络图Returns-------None.'''# ------------------ all inputs ------------------------self.s = tf.placeholder(tf.float32, [None, self.n_features], name='s')  # input Stateself.s_ = tf.placeholder(tf.float32, [None, self.n_features], name='s_')  # input Next Stateself.r = tf.placeholder(tf.float32, [None, ], name='r')  # input Rewardself.a = tf.placeholder(tf.int32, [None, ], name='a')  # input Action# 满足正态分布w_initializer, b_initializer = tf.random_normal_initializer(0., 0.3), tf.constant_initializer(0.1)# ------------------ build evaluate_net ------------------with tf.variable_scope('eval_net'):# tf.layers.dense()：添加一个全连接层# tf.nn.relu6：计算校正线性6：min(max(features, 0), 6)e1 = tf.layers.dense(self.s, 100, tf.nn.relu6, kernel_initializer=w_initializer,bias_initializer=b_initializer, name='e1')# e2 = tf.layers.dense(e1, 48, tf.nn.relu6, kernel_initializer=w_initializer,#                      bias_initializer=b_initializer, name='e2')e3 = tf.layers.dense(e1, 20, tf.nn.relu, kernel_initializer=w_initializer,bias_initializer=b_initializer, name='e3')self.q_eval = tf.layers.dense(e3, self.n_actions, tf.nn.softmax, kernel_initializer=w_initializer,bias_initializer=b_initializer, name='q')# ------------------ build target_net ------------------with tf.variable_scope('target_net'):t1 = tf.layers.dense(self.s_, 100, tf.nn.relu6, kernel_initializer=w_initializer,bias_initializer=b_initializer, name='t1')# t2 = tf.layers.dense(t1, 48, tf.nn.relu6, kernel_initializer=w_initializer,#                      bias_initializer=b_initializer, name='t2')t3 = tf.layers.dense(t1, 20, tf.nn.relu, kernel_initializer=w_initializer,bias_initializer=b_initializer, name='t3')self.q_next = tf.layers.dense(t3, self.n_actions, tf.nn.softmax, kernel_initializer=w_initializer,bias_initializer=b_initializer, name='t4')with tf.variable_scope('q_target'):q_target = self.r + self.gamma * tf.reduce_max(self.q_next, axis=1, name='Qmax_s_')  # shape=(None, )self.q_target = tf.stop_gradient(q_target)with tf.variable_scope('q_eval'):a_indices = tf.stack([tf.range(tf.shape(self.a)[0], dtype=tf.int32), self.a], axis=1)self.q_eval_wrt_a = tf.gather_nd(params=self.q_eval, indices=a_indices)  # shape=(None, )with tf.variable_scope('loss'):self.loss = tf.reduce_mean(tf.squared_difference(self.q_target, self.q_eval_wrt_a, name='TD_error'))with tf.variable_scope('train'):# self._train_op = tf.train.RMSPropOptimizer(self.lr).minimize(self.loss)self._train_op = tf.train.AdamOptimizer(self.lr).minimize(self.loss)def store_transition(self, s, a, r, s_):'''存储记忆元组Parameters----------s : 当前状态a : 动作r : 奖惩s_ : 下一状态Returns-------None.'''if not hasattr(self, 'memory_counter'):self.memory_counter = 0# np.hstack将参数元组的元素数组按水平方向进行叠加transition = np.hstack((s, a, [r], s_))# replace the old memory with new memoryindex = self.memory_counter % self.memory_sizeself.memory[index, :] = transitionself.memory_counter += 1def choose_action(self, observation):'''根据当前状态选择动作Parameters----------observation : 对环境的观测Returns-------action : 做出的动作'''# to have batch dimension when feed into tf placeholder# 环境状态observation = observation[np.newaxis, :]# e贪心算法if np.random.uniform() < self.epsilon:# forward feed the observation and get q value for every actionsactions_value = self.sess.run(self.q_eval, feed_dict={self.s: observation})action = np.argmax(actions_value)else:action = np.random.randint(0, self.n_actions)return actiondef learn(self):'''神经网络训练Returns-------None.'''# check to replace target parametersif self.learn_step_counter % self.replace_target_iter == 0:self.sess.run(self.target_replace_op)print('\ntarget_params_replaced\n')# sample batch memory from all memory# 从所有记忆中选取小批量记忆if self.memory_counter > self.memory_size:sample_index = np.random.choice(self.memory_size, size=self.batch_size)else:sample_index = np.random.choice(self.memory_counter, size=self.batch_size)batch_memory = self.memory[sample_index, :]_, cost = self.sess.run([self._train_op, self.loss],feed_dict={self.s: batch_memory[:, :self.n_features],self.a: batch_memory[:, self.n_features],self.r: batch_memory[:, self.n_features + 1],self.s_: batch_memory[:, -self.n_features:],})self.cost_his.append(cost)# increasing epsilon# self.epsilon = self.epsilon + self.epsilon_increment if self.epsilon < self.epsilon_max else self.epsilon_maxself.learn_step_counter += 1