-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathfunction_wmmse_powercontrol.py
172 lines (157 loc) · 7.01 KB
/
function_wmmse_powercontrol.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
# ###############################################
# This file includes functions to perform the WMMSE algorithm [2].
# Codes have been tested successfully on Python 3.6.0 with Numpy 1.12.0 support.
#
# References:
# [1] Haoran Sun, Xiangyi Chen, Qingjiang Shi, Mingyi Hong, Xiao Fu, and Nikos D. Sidiropoulos.
# "Learning to optimize: Training deep neural networks for wireless resource management."
# in proceedings of IEEE 18th International Workshop on Signal Processing Advances in Wireless Communications (SPAWC), 2017.
#
# [2] Haoran Sun, Xiangyi Chen, Qingjiang Shi, Mingyi Hong, Xiao Fu, and Nikos D. Sidiropoulos,
# "Learning to Optimize: Training Deep Neural Networks for Interference Management,"
# in IEEE Transactions on Signal Processing, vol. 66, no. 20, pp. 5438-5453, 15 Oct.15, 2018.
#
# [3] Qingjiang Shi, Meisam Razaviyayn, Zhi-Quan Luo, and Chen He.
# "An iteratively weighted MMSE approach to distributed sum-utility maximization for a MIMO interfering broadcast channel."
# IEEE Transactions on Signal Processing 59, no. 9 (2011): 4331-4340.
#
# version 1.0 -- February 2017. Written by Haoran Sun (hrsun AT iastate.edu)
# ###############################################
import numpy as np
import math
import time
import scipy.io as sio
import matplotlib.pyplot as plt
# Functions for objective (sum-rate) calculation
def obj_IA_sum_rate(H, p, var_noise, K):
y = 0.0
for i in range(K):
s = var_noise
for j in range(K):
if j!=i:
s = s+H[i,j]**2*p[j]
y = y+math.log2(1+H[i,i]**2*p[i]/s)
return y
def obj_IA_sum_log_rate(H, p, var_noise, K):
y = 0.0
for i in range(K):
s = var_noise
for j in range(K):
if j!=i:
s = s+H[i,j]**2*p[j]
y = y+math.log10(math.log2(1+H[i,i]**2*p[i]/s))
return y
# Functions for WMMSE algorithm
def WMMSE_sum_rate(p_int, H, Pmax, var_noise):
K = np.size(p_int)
vnew = 0
b = np.sqrt(p_int)
f = np.zeros(K)
w = np.zeros(K)
for i in range(K):
f[i] = H[i, i] * b[i] / (np.square(H[i, :]) @ np.square(b) + var_noise)
w[i] = 1 / (1 - f[i] * b[i] * H[i, i])
vnew = vnew + math.log2(w[i])
VV = np.zeros(10000)
for iter_ in range(10000):
vold = vnew
for i in range(K):
btmp = w[i] * f[i] * H[i, i] / sum(w * np.square(f) * np.square(H[:, i]))
b[i] = min(btmp, np.sqrt(Pmax)) + max(btmp, 0) - btmp
vnew = 0
for i in range(K):
f[i] = H[i, i] * b[i] / ((np.square(H[i, :])) @ (np.square(b)) + var_noise)
w[i] = 1 / (1 - f[i] * b[i] * H[i, i])
vnew = vnew + math.log2(w[i])
VV[iter_] = vnew
if vnew - vold <= 1e-5:
break
p_opt = np.square(b)
return p_opt, iter_
# Functions for performance evaluation
def perf_eval(H, Py_p, NN_p, K, var_noise=1):
num_sample = H.shape[2]
pyrate = np.zeros(num_sample)
nnrate = np.zeros(num_sample)
mprate = np.zeros(num_sample)
rdrate = np.zeros(num_sample)
for i in range(num_sample):
pyrate[i] = obj_IA_sum_rate(H[:, :, i], Py_p[:, i], var_noise, K)
nnrate[i] = obj_IA_sum_rate(H[:, :, i], NN_p[i, :], var_noise, K)
mprate[i] = obj_IA_sum_rate(H[:, :, i], np.ones(K), var_noise, K)
rdrate[i] = obj_IA_sum_rate(H[:, :, i], np.random.rand(K,1), var_noise, K)
print('Sum-rate: WMMSE: %0.3f, DNN: %0.3f, Max Power: %0.3f, Random Power: %0.3f'%(sum(pyrate)/num_sample, sum(nnrate)/num_sample, sum(mprate)/num_sample, sum(rdrate)/num_sample))
print('Ratio: DNN: %0.3f%%, Max Power: %0.3f%%, Random Power: %0.3f%%\n' % (sum(nnrate) / sum(pyrate)* 100, sum(mprate) / sum(pyrate) * 100, sum(rdrate) / sum(pyrate) * 100))
plt.figure('%d'%K)
plt.style.use('seaborn-deep')
data = np.vstack([pyrate, nnrate]).T
bins = np.linspace(0, max(pyrate), 50)
plt.hist(data, bins, alpha=0.7, label=['WMMSE', 'DNN'])
plt.legend(loc='upper right')
plt.xlim([0, 8])
plt.xlabel('sum-rate')
plt.ylabel('number of samples')
plt.savefig('Histogram_%d.eps'%K, format='eps', dpi=1000)
plt.show()
return 0
# Functions for data generation, Gaussian IC case
def generate_Gaussian(K, num_H, Pmax=1, Pmin=0, seed=2017):
print('Generate Data ... (seed = %d)' % seed)
np.random.seed(seed)
Pini = Pmax*np.ones(K)
var_noise = 1
X=np.zeros((K**2,num_H))
Y=np.zeros((K,num_H))
total_time = 0.0
for loop in range(num_H):
CH = 1/np.sqrt(2)*(np.random.randn(K,K)+1j*np.random.randn(K,K))
H=abs(CH)
X[:,loop] = np.reshape(H, (K**2,), order="F")
H=np.reshape(X[:,loop], (K,K), order="F")
mid_time = time.time()
Y[:,loop] = WMMSE_sum_rate(Pini, H, Pmax, var_noise)
total_time = total_time + time.time() - mid_time
# print("wmmse time: %0.2f s" % total_time)
return X, Y, total_time
# Functions for data generation, Gaussian IC half user case
def generate_Gaussian_half(K, num_H, Pmax=1, Pmin=0, seed=2017):
print('Generate Testing Data ... (seed = %d)' % seed)
np.random.seed(seed)
Pini = Pmax * np.ones(K)
var_noise = 1
X = np.zeros((K ** 2 * 4, num_H))
Y = np.zeros((K * 2, num_H))
total_time = 0.0
for loop in range(num_H):
CH = 1 / np.sqrt(2) * (np.random.randn(K, K) + 1j * np.random.randn(K, K))
H = abs(CH)
print('H is',H)
print('Pini is', Pini)
mid_time = time.time()
Y[0: K, loop] = WMMSE_sum_rate(Pini, H, Pmax, var_noise)
total_time = total_time + time.time() - mid_time
OH = np.zeros((K * 2, K * 2))
OH[0: K, 0:K] = H
X[:, loop] = np.reshape(OH, (4 * K ** 2,), order="F")
# print("wmmse time: %0.2f s" % total_time)
return X, Y, total_time ##X is the channel matrix vectorized to form a column for the NN,
## If you set K = 5 (number of users), X will have a length of 5x5
## Y is vector of the power allocated to the 5 users which will serve
## as a label for the NN. num_H is the number of dataset like snapshots
# Functions for data generation, IMAC case
def generate_IMAC(num_BS, num_User, num_H, Pmax=1, var_noise = 1):
# Load Channel Data
CH = sio.loadmat('IMAC_%d_%d_%d' % (num_BS, num_User, num_H))['X']
Temp = np.reshape(CH, (num_BS, num_User * num_BS, num_H), order="F")
H = np.zeros((num_User * num_BS, num_User * num_BS, num_H))
for iter in range(num_BS):
H[iter * num_User:(iter + 1) * num_User, :, :] = Temp[iter, :, :]
# Compute WMMSE output
Y = np.zeros((num_User * num_BS, num_H))
Pini = Pmax * np.ones(num_User * num_BS)
start_time = time.time()
for loop in range(num_H):
Y[:, loop] = WMMSE_sum_rate(Pini, H[:, :, loop], Pmax, var_noise)
wmmsetime=(time.time() - start_time)
# print("wmmse time: %0.2f s" % wmmsetime)
return CH, Y, wmmsetime, H