目录
MEC网络，以时间间隔来划分整个网络的运行过程，包括三种行为，计算卸载(WD->AP)和充能过程(AP-WD)，行为在一个时间间隔内进行，Download行为是返回卸载的结果 ，这是一个时分复用的网络，所有的行为都在一个频带内进行，因此对于一个ap下的所有节点，必须规划整个的过程。
整个问题为最大化计算速率，肢解成两个部分，包括问题1，计算卸载策略，问题二，资源分配。问题一用深度强化学习，问题二用一维双截面搜索。
源码描述之optimization.py
用于解决资源分配问题，这是一个凸优化问题
此时卸载策略已经确定，1为卸载，0为local计算
主要是关于这个文件的注释，剩下两个文件一个memoryDNN是训练需要的神经网络。另一个是训练卸载策略的主文件，这个文件才是问题解决的核心。
# -*- coding: utf-8 -*-
"""
Created on Tue Jan 9 10:45:26 2018
@author: Administrator
"""
import numpy as np
from scipy import optimize
from scipy.special import lambertw
import scipy.io as sio # import scipy.io for .mat file I/
import time
def plot_gain( gain_his):
import matplotlib.pyplot as plt
import pandas as pd
import matplotlib as mpl
gain_array = np.asarray(gain_his)
df = pd.DataFrame(gain_his)
mpl.style.use('seaborn')
fig, ax = plt.subplots(figsize=(15,8))
rolling_intv = 20
plt.plot(np.arange(len(gain_array))+1, df.rolling(rolling_intv, min_periods=1).mean(), 'b')
plt.fill_between(np.arange(len(gain_array))+1, df.rolling(rolling_intv, min_periods=1).min()[0], df.rolling(rolling_intv, min_periods=1).max()[0], color = 'b', alpha = 0.2)
plt.ylabel('Gain ratio')
plt.xlabel('learning steps')
plt.show()
#线性搜索函数，需要两个输入，一个是信道矩阵，另一个是卸载动作，关于每个无线设备，除了W参数，其他的默认是同质的（我对代码的理解）。
#这里的计算过程来自于上面问题求解的公式
def bisection(h, M, weights=[]):
# the bisection algorithm proposed by Suzhi BI
# average time to find the optimal: 0.012535839796066284 s
# parameters and equations
o=100
#是处理一位原始数据所需的周期数
p=3
#AP功率
u=0.7
#能量收集效率
eta1=((u*p)**(1.0/3))/o
#一个固定参数，不用理解
ki=10**-26
#能量效率系数，跟计算机能耗相关的
eta2=u*p/10**-10
B=2*10**6
#带宽
Vu=1.1
#Vu是原始数据（卸载数据）经过加密后的数据（原有数据在传输时进行额外添加）大小比原始数据的倍数
epsilon=B/(Vu*np.log(2))
x = [] # a =x[0], and tau_j = a[1:],时间分配矩阵
M0=np.where(M==0)[0]
#为零或者1的索引
M1=np.where(M==1)[0]
hi=np.array([h[i] for i in M0])
#1或者0索引对应的信道
hj=np.array([h[i] for i in M1])
if len(weights) == 0:
# default weights [1, 1.5, 1, 1.5, 1, 1.5, ...]
weights = [1.5 if i%2==1 else 1 for i in range(len(M))]
wi=np.array([weights[M0[i]] for i in range(len(M0))])
wj=np.array([weights[M1[i]] for i in range(len(M1))])
#w是每个WD的权重
def sum_rate(x):
sum1=sum(wi*eta1*(hi/ki)**(1.0/3)*x[0]**(1.0/3))
sum2=0
for i in range(len(M1)):
sum2+=wj[i]*epsilon*x[i+1]*np.log(1+eta2*hj[i]**2*x[0]/x[i+1])
return sum1+sum2
def phi(v, j):
return 1/(-1-1/(lambertw(-1/(np.exp( 1 + v/wj[j]/epsilon))).real))
def p1(v):
p1 = 0
for j in range(len(M1)):
p1 += hj[j]**2 * phi(v, j)
return 1/(1 + p1 * eta2)
def Q(v):
sum1 = sum(wi*eta1*(hi/ki)**(1.0/3))*p1(v)**(-2/3)/3
#local计算
sum2 = 0
for j in range(len(M1)):
sum2 += wj[j]*hj[j]**2/(1 + 1/phi(v,j))
#计算卸载
return sum1 + sum2*epsilon*eta2 - v
#这里来自公式19
def tau(v, j):
return eta2*hj[j]**2*p1(v)*phi(v,j)
# bisection starts here
delta = 0.005
UB = 999999999
LB = 0
while UB - LB > delta:
v = (float(UB) + LB)/2
if Q(v) > 0:
LB = v
else:
UB = v
x.append(p1(v))
#计算a
for j in range(len(M1)):
x.append(tau(v, j))
#计算每个WD的t
return sum_rate(x), x[0], x[1:]
#CD法的代码
def cd_method(h):
N = len(h)
M0 = np.random.randint(2,size = N)
gain0,a,Tj= bisection(h,M0)
g_list = []
M_list = []
while True:
for j in range(0,N):
M = np.copy(M0)
M[j] = (M[j]+1)%2
gain,a,Tj= bisection(h,M)
g_list.append(gain)
M_list.append(M)
g_max = max(g_list)
if g_max > gain0:
gain0 = g_max
M0 = M_list[g_list.index(g_max)]
else:
break
return gain0, M0
if __name__ == "__main__":
h=np.array([6.06020304235508*10**-6,1.10331933767028*10**-5,1.00213540309998*10**-7,1.21610610942759*10**-6,1.96138838395145*10**-6,1.71456339592966*10**-6,5.24563569673585*10**-6,5.89530717142197*10**-7,4.07769429231962*10**-6,2.88333185798682*10**-6])
M=np.array([1,0,0,0,1,0,0,0,0,0])
# h=np.array([1.00213540309998*10**-7,1.10331933767028*10**-5,6.06020304235508*10**-6,1.21610610942759*10**-6,1.96138838395145*10**-6,1.71456339592966*10**-6,5.24563569673585*10**-6,5.89530717142197*10**-7,4.07769429231962*10**-6,2.88333185798682*10**-6])
# M=np.array([0,0,1,0,1,0,0,0,0,0])
# h = np.array([4.6368924987170947*10**-7, 1.3479411763648968*10**-7, 7.174945246007612*10**-6, 2.5590719803595445*10**-7, 3.3189928740379023*10**-6, 1.2109071327755575*10**-5, 2.394278475886022*10**-6, 2.179121774067472*10**-6, 5.5213902658478367*10**-8, 2.168778154948169*10**-7, 2.053227965874453*10**-6, 7.002952297466865*10**-8, 7.594077851181444*10**-8, 7.904048961975136*10**-7, 8.867218892023474*10**-7, 5.886007653360979*10**-6, 2.3470565740563855*10**-6, 1.387049627074303*10**-7, 3.359475870531776*10**-7, 2.633733784949562*10**-7, 2.189895264149453*10**-6, 1.129177795302099*10**-5, 1.1760290137191366*10**-6, 1.6588656719735275*10**-7, 1.383637788476638*10**-6, 1.4485928387351664*10**-6, 1.4262265958416598*10**-6, 1.1779725004265418*10**-6, 7.738218993031842*10**-7, 4.763534225174186*10**-6])
# M =np.array( [0, 0, 1, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 1,])
# time the average speed of bisection algorithm
# repeat = 1
# M =np.random.randint(2, size=(repeat,len(h)))
# start_time=time.time()
# for i in range(repeat):
# gain,a,Tj= bisection(h,M[i,:])
# total_time=time.time()-start_time
# print('time_cost:%s'%(total_time/repeat))
gain,a,Tj= bisection(h,M)
print('y:%s'%gain)
print('a:%s'%a)
print('Tj:%s'%Tj)
# test CD method. Given h, generate the max mode
#就是Cd法，可以见关于cd法的介绍
#显然在user少的时候CD法很快，但随着用户数量增加，时间开销会很大
gain0, M0 = cd_method(h)
print('max y:%s'%gain0)
print(M0)
# test all data
K = [10, 20, 30] # number of users
N = 1000 # number of channel
#对于数据集合中的数据使用CD法来优化
for k in K:
# Load data
channel = sio.loadmat('./data/data_%d' %int(k))['input_h']
gain = sio.loadmat('./data/data_%d' %int(k))['output_obj']
start_time=time.time()
gain_his = []
gain_his_ratio = []
mode_his = []
for i in range(N):
if i % (N//10) == 0:
print("%0.1f"%(i/N))
i_idx = i
h = channel[i_idx,:]
# the CD method
gain0, M0 = cd_method(h)
# memorize the largest reward
gain_his.append(gain0)
gain_his_ratio.append(gain_his[-1] / gain[i_idx][0])
mode_his.append(M0)
total_time=time.time()-start_time
print('time_cost:%s'%total_time)
print('average time per channel:%s'%(total_time/N))
plot_gain(gain_his_ratio)
print("gain/max ratio: ", sum(gain_his_ratio)/N)
这串代码涉及的公式包括CD法来自于论文文献参考的第五篇，但我下面的概述已经足以解释了，如果想看详细的证明过程可以去看那篇论文。
公式描述
本文主要描述一个类似DROO的场景，主要在于提出一种一维的优化搜索算法。介绍本文的建模：
无线信道的能量收集：
p是发射功率，u是接收比例，h是信道增益，T是一个时间间隙，a是所占时间间隙的比例
local计算：
对于local来说，计算和能量接收时可以并行的，因为两个存在于不同的电路中，local可以在整个T的时间内进行计算。
能量上限和周期数的关系，
这里不考虑异构性的影响，对于某个计算任务来说，在不同的WD上计算开销都是相同的，
是处理一位原始数据所需的周期数，这里不对能量变化做考虑，也就是说WD总是可以维持一个功率来一直消耗完所有能量。
f代表处理器的每秒周期数，t是WD的计算时间，r是local的计算速度，bit/s
这里我大胆推测一下，显然保持低的每秒时钟周期数可以使得最后的ft最大，也就是单个MD在T时间内计算的数据更多，也就是说要使得另一个变量t达到最大，也就是T，因此，可以得到的结果是。
卸载到边：
必须在能量收集阶段之后，
卸载时的信道容量
这里的b代表了可以卸载多少原始数据到AP上，信道每秒的
传输量，tT是传输时间，Vu是原始数据经过加密后的数据大小比原始数据的倍数，因此这里b代表在这种情况下，时间内能卸载原始数据的数据量。
这是卸载原始数据b对应的边计算并返回的时间，前面是计算时间，后面是返回时间。rd<<1,代表返回数据和原始数据b的比例。因此可以说li在总体时间计算时可以忽略，而只计算tT和aT的时间，这是一个时隙内的主要时间开销。
这是最后mode1下WD的传输速率，下面用的是WD的传输速率来代表WD的计算速率，耗能采取尽量全部耗尽的方式（在分配的时间内用最大的功率来传输，不设置上限，因为考虑接收到的能量是很少的），而WD的作用就是尽力传输任务到边节点。
两个计算速率都是相对于T而言的，表示时间间隔为T的时候WD的计算卸载速率，其实就是在时间间隔内WD可以计算的数据的总量。
所有的节点的综合的计算速率：
w是每个WD的参数。
local：
这里有，本来还自己算来着，
，这个参数被固定了，因此问题的式的c和d两个条件也不需要了
。
问题一的简化版本，根据local的最大化的推理式得到，可以看到如果卸载策略给定之前，需要确定一个卸载策略，然后在这个策略下确定一个最优的资源分配，显然这是一个非凸的问题，而如果策略给定之后，这就是一个凸优化的问题，而目标是找到a和t的分配比例。
后面的主要工作是在给定卸载策略后，传输时间的分配。
问题的拉格朗日数乘法：
对偶问题是
而在最优解处，
的等号是一定成立的。
三个参数之间的关系，证明没看，知道就行。在推导处v大于0，因为分母被证明是一定大于0的，而且w(x)在（1/e，0）是增函数，因此，可以判断更长的t一定是跟大的h和w参数相关，一个是无线信道强，一个是权重大。
这是t和a关系的函数，
带入t+a =1，两者和为一的等式：
得到a和v的关系式，
得到v的关系式，这个式来自对于a的偏导，带入a = p(v)，因此得到求v值的目标函数式，使用二分搜索来求这个唯一的v值，而后根据v来求a，然后求t。
求v的二分搜索算法，然后计算a和t
上述式第一种方法，解耦整个优化问题，分成两个部分优化，使用二分搜索算法求解拉格朗日方程。这种方式需要先假定计算模式选择，然后优化传输时间分配，当WD设备很多时，这种方式使得算法的复杂度会很大。
CD法
然后使用坐标下降法来优化模式选择，这是一个nphard问题
一次优化一个m
算法，大致意思就是每轮都更换每个wd的模式，然后挑选max的奖励，
结束是更换任意一个WD的模式都没有一个正向的奖励，可以说是一个极大值解，但不一定是最优解。
这里有一个同质的示例，当每个WD相同时，因素只剩信道增益，而信道增益大的WD将更多会卸载到边上。
下面用另一个方法：
交替方向乘子法ADMM
用于联合优化计算模式选择和传输时间发分配，这里把硬组合优化问题分解成N个并行较小的整数规划问题，每个WD一个，但WD之间不是独立的，参数是耦合的；引入二元决策变量m把优化问题的两个部分合二为一。
算法流程
大致的意思是相同的，不过在这里把m和a和t放在同一个级别上，并且对他们进行交替的求最大值然后更新参数。
比单纯的test每个m的计算要少很多。
DROO的网络模型
包括两部分，一部分是深度网络，用来生成卸载决策，一部分则根据决策求解两个值，完成之后根据奖励对策略进行更新。
略谈边缘计算卸载模型
卸载模式
二进制卸载和部分卸载，二进制卸载用于不可分割的简单任务，对于一个任务来说，只有两种选择，卸载or不卸载，而部分卸载则有利于一些由多个并行段组成的复杂任务，部分卸载又分为平行和顺序结构，对于平行结构来说，任务的组件之间是独立的，对于顺序结构来说，任务之间的执行具有前后关系，如果是全粒度划分，任务可以以任意的比例被分割成不同的部分，一个任务可以以任意的比例分割成两个任务，具有非全粒度任务分割的部分卸载可以看作二进制卸载，因此整数规划在卸载建模过程中经常被用到。
信道模型
FDMA
TDMA
CDMA
SDMA
OFDMA
计算模型
两点： 时间和能量，
本地计算的时延和能量开销，数据传输的能量消耗和延迟，边缘计算的能量开销和延迟
任务的抽象方式：消耗周期数大小，任务数据大小+数据和周期消耗的比例
能量收集模型
确定性模型，随机模型
单位时间内提供的能量是已知的，例如定向的充电
或者能量得到是随机的，基于 某种模型
关于np问题的解释
np问题：多项式时间很难找到解，但解可以在多项式时间内验证。
npc问题：所有np问题都可以规划到一个npc问题。
np-hard，验证解也无法在多项式时间内做到。
了解这个问题后关于强化学习的思考
应该考虑的是把原本就适合使用强化学习的问题进而再使用dl来进行结合。
已有的一些成果的原因往往是因为针对问题的性质，适宜强化学习。
DQN和AlphaGo系列工作给人留下深刻印象，但是这两种任务本质上其实相对“简单”。因为这些任务的环境是确定的、静态的，状态主要是离散的、静态的、完全可观察的，反馈是确定的，代理也是单一的。目前DRL在解决部分可见状态任务（如StarCraft），状态连续的任务（如机械控制任务），动态反馈任务和多代理任务中还没取得令人惊叹的突破。