-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy pathErlangB_Simulation.py
184 lines (158 loc) · 6.07 KB
/
ErlangB_Simulation.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
173
174
175
176
177
178
179
180
181
182
183
184
from scipy.stats import lognorm
from scipy.stats import expon
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
import random
import math
#max number of calls handled
n = 41
#values for random call generation
maxValue = 3600
skewness = 0.5
average_no_calls=99
std_deviation = 27
graph_x=[]
graph_y=[]
graph_x2=[]
graph_y2=[]
#generate random varibales for simulation
def create_random_variables():
#generate random # of calls from normal distribution
no_calls = int(random.normalvariate(average_no_calls,std_deviation))
if no_calls <= 0:
#if no_calls invalid set to standard avergae # calls not in busy hour
no_calls = 67
#get random length of calls
#random_call_length2 = lognorm.rvs(s =skewness,loc=average_no_calls, size=no_calls)
random_call_length2 = np.random.random_sample(size = no_calls)
c = 0
while c < len(random_call_length2):
x = 1 - random_call_length2[c]
random_call_length2[c] = math.log(x)/ -(1/900)
c+=1
#checks to make sure call length is ok
if (len(random_call_length2) == 0):
random_call_length2 = [0.1]
for call in random_call_length2:
if call < 0 or np.isnan(call):
random_call_length2[call] = 0.1
# np.seterr(all='raise')
# #normalise data between 0-3600 seconds
# try:
# random_call_length2 = random_call_length2 - min(random_call_length2)
# random_call_length2 = random_call_length2 / max(random_call_length2)
# random_call_length2 = random_call_length2 * maxValue
# except FloatingPointError:
# print("Invalid value caught and programme continues")
#put random call lengths into a list parsed into ints
random_call_lengths_int = []
for value in random_call_length2:
if math.isnan(value):
random_call_lengths_int.append(900)
continue
random_call_lengths_int.append(int(value))
# plt.hist(random_call_length2,30,density=True, color = 'red', alpha=0.1)
# plt.xlabel("Duration in Seconds (s)")
# plt.ylabel("Probability of occurance")
# plt.show()
p = 0
time_between_calls = int(maxValue/no_calls)
time=0
call_dict = {}
#assign call length to call start time
while p < no_calls:
time += time_between_calls
call_dict.update({time:random_call_lengths_int[p]})
p+=1
return call_dict
#carry out Monte Carlo Simulation
def simulate_calls(call_dict):
global graph_x
global graph_y
time = 0
simultaneous_calls=0
dropped_calls={}
skip_because_dropped=[]
while time < maxValue:
for call in call_dict.items():
#If time = call start time, add it to simultaneous calls count, given there is room
#otherwise drop call and add it to dropped_calls list
if call[0] == time:
if simultaneous_calls < n:
simultaneous_calls +=1
else:
dropped_calls.update({call[0]:call[1]})
skip_because_dropped.append(call[0])
#If call end time = time, take it off simultaneous call count
#Check to make sure that call has not been dropped
if (call[1] + call[0]) == time:
if call[0] in skip_because_dropped:
continue
simultaneous_calls -=1
time +=1
#get params to calculate GOS manually
avg_call_duration = sum(call_dict.values())/len(call_dict.values())
avg_offered_traffic = (avg_call_duration/maxValue)*len(call_dict.values())
avg_GOS_from_dropped = 100*len(dropped_calls)/len(call_dict)
#add Ao and GOS to the list for graphing
graph_x.append(avg_offered_traffic)
graph_y.append(avg_GOS_from_dropped)
return(avg_call_duration,len(call_dict),avg_offered_traffic,avg_GOS_from_dropped,len(dropped_calls))
def calculate_GOS_erlangb(bottom_range,top_range):
global graph_x2
global graph_y2
all_data= []
iterator = bottom_range
while iterator <= top_range:
Ao = iterator
i = 1
n_factorial = 1
factorial_list=[]
numerator = 0
#calculate factorial
while i <= n:
n_factorial= n_factorial * i
factorial_list.append(n_factorial)
i +=1
#calculate numerator
numerator = (Ao**n)/factorial_list[n-1]
denominator=0
j = 1
#calculate denominator
while j <= len(factorial_list):
denominator +=(Ao**j)/factorial_list[j-1]
j+=1
denominator +=1
#calculate GOS
E1 = numerator/denominator
iterator +=1
all_data.append([Ao,E1*100])
graph_x2.append(Ao)
graph_y2.append(E1*100)
return all_data
if __name__ == "__main__":
simulation_list = []
print("Please wait a moment while the programme executes...")
k=0
while k < 1000:
call_dictionary = create_random_variables()
call_dur,calls,offered_traffic,GOS,dropped_calls =simulate_calls(call_dictionary)
simulation_list.append([call_dur,calls,offered_traffic,GOS,dropped_calls])
k+=1
print("\nResults from Monte deCarlo simulation with Erlang B.\nOffered traffic varied with random number of calls and varied call length.\nConstant channels = 41")
results_df = pd.DataFrame.from_records(simulation_list, columns=['Avg Call Duration',
'Avg No Calls',
'Avg Offered Traffic',
'Avg GOS',
'Avg No Dropped Calls'])
print(results_df.describe())
results_df.describe().style.format('{:,}')
xs, ys = zip(*sorted(zip(graph_x, graph_y)))
calculate_GOS_erlangb(results_df['Avg Offered Traffic'].min(),results_df['Avg Offered Traffic'].max())
plt.plot(graph_x2,graph_y2)
plt.plot(xs,ys)
plt.xlabel("Ao")
plt.ylabel("GOS")
plt.title("Erlang B formula Vs Simulation with Decaying Exponential")
plt.show()