1
+ %% This script simulates the Lyapunov Optimization-based Dynamic Computation Offloading (LODCO) algorithm.
2
+ % author: Hailiang Zhao
3
+ clc , clear
4
+ opt = optimset(' Display' , ' none' );
5
+
6
+ %% basic parameter settings
7
+ k = 1e- 28 ; % effective switched capacitance (a constant decided by the chip architecture)
8
+ tau = 0.002 ; % the length of time slot (in second)
9
+ phi = 0.002 ; % the cost of task dropping (in second)
10
+ omega = 1e6 ; % the bandwidth of MEC server (in Hz)
11
+ sigma = 1e- 13 ; % the noise power of the receiver (in W)
12
+ p_tx_max = 1 ; % the maximum transmit power of mobile device (in W)
13
+ f_max = 1.5e9 ; % the maximum CPU-cycle frequency of mobile device (in Hz)
14
+ E_max = 0.002 ; % the maximum amout of battery output energy (in J)
15
+ L = 1000 ; % the input size of the computation task (in bit)
16
+ X = 737.5 ; % the number of CPU cycles needed on processing one bit of task
17
+ W = L * X ; % the number of CPU cycles needed on processing one task
18
+ E_H_max = 48e-6 ; % the upper bound of the energy arrive at the mobile device (in J)
19
+ p_H = E_H_max / (2 * tau ); % the average Energy Harvesting (EH) power (in W)
20
+ g0 = power(10 , - 4 ); % the path-loss constant
21
+
22
+ %% parameter control
23
+ T = 10000 ; % the number of time slot (a.k.a. the size of the time horizon)
24
+ tau_d = 0.002 ; % execution deadline (in second)
25
+ d = 50 ; % the distance between the mobile device and the MEC server (in meter)
26
+ E_min = 0.02e-3 ; % the minimum amout of battery output energy (in J)
27
+ V = 1e-5 ; % the weight of penalty (the control parameter introduced by Lyapunov Optimization)
28
+ rho = 0.6 ; % the probability that the computation task is requested
29
+
30
+ % the lower bound of perturbation parameter
31
+ E_max_hat = min(max(k * W * (f_max )^2 , p_tx_max * tau ), E_max );
32
+ theta = E_max_hat + V * phi / E_min ;
33
+
34
+ %% allocate storage for valuable results
35
+ B = zeros(T , 1 ); % the battery energy level (in J)
36
+ B_hat = zeros(T , 1 ); % the virtual battery energy level ($B_hat = B - theta$)
37
+ e = zeros(T , 1 ); % the amout of the harvested and stored energy (in J)
38
+ chosen_mode = zeros(T , 1 );% {1: local, 2: remote, 3: drop, 4: no task request}
39
+ f = zeros(T , 1 ); % the CPU-cycle frequency of local execution (in Hz)
40
+ p = zeros(T , 1 ); % the transmit power of computation offloading (in W)
41
+ cost = zeros(T , 3 ); % execution delay for mobile execution, MEC server execution and final choice, respectively (in second)
42
+ E = zeros(T , 3 ); % energy consumption for mobile execution, MEC server execution and final choice, respectively (in J)
43
+
44
+ %% simulation begin
45
+ t = 1 ;
46
+ while t <= T
47
+ disp([' ===> Time slot #' , num2str(t ), ' <===' ])
48
+
49
+ %% initialization
50
+ % generate the task request
51
+ zeta = binornd(1 , rho );
52
+ % generate the virtual battery energy level
53
+ B_hat(t ) = B(t ) - theta ;
54
+
55
+ %% step 1: get the optimal energy harvesting no matter whether task is requested
56
+ E_H_t = unifrnd(0 , E_H_max );
57
+ if B_hat(t ) <= 0
58
+ e(t ) = E_H_t ;
59
+ end
60
+
61
+ %% step 2: get the optimal computation offloading strategy (I_m, I_s, I_d, f(t), p(t))
62
+ if zeta == 0
63
+ % chosen mode has to be 4
64
+ disp(' no task request generated!' )
65
+ chosen_mode(t ) = 4 ;
66
+ else
67
+ % chosen_mode is chosen from {1, 2, 3}
68
+ disp(' task request generated!' );
69
+ % task request exists, generate the channel power gain
70
+ h = exprnd(g0 / power(d , 4 ));
71
+
72
+ %% step 2.1: solve the optimization problem $\mathcal{P}_{ME}$ (f(t) > 0)
73
+ % calculate f_L and f_U
74
+ f_L = max(sqrt(E_min / (k * W )), W / tau_d );
75
+ f_U = min(sqrt(E_max / (k * W )), f_max );
76
+ if f_L <= f_U
77
+ % the sub-problem is feasible
78
+ disp(' mobile execution ($\mathcal{P}_{ME}$) is feasible!' )
79
+
80
+ if B_hat(t ) < 0
81
+ f_0 = (V / (-2 * B_hat(t ) * k ))^(1 / 3 );
82
+ else
83
+ % complex number may exist, which will lead to error
84
+ f_0 = -(V / (2 * B_hat(t ) * k ))^(1 / 3 );
85
+ end
86
+
87
+ if (f_0 > f_U && B_hat(t ) < 0 ) || (B_hat(t ) >= 0 )
88
+ f(t ) = f_U ;
89
+ elseif f_0 >= f_L && f_0 <= f_U && B_hat(t ) < 0
90
+ f(t ) = f_0 ;
91
+ elseif f_0 < f_L && B_hat(t ) < 0
92
+ f(t ) = f_L ;
93
+ end
94
+ % check whether f(t) is zero
95
+ if f(t ) == 0
96
+ disp(' Something wrong! f is 0!' )
97
+ end
98
+
99
+ % calculate the delay of mobile execution
100
+ cost(t , 1 ) = W / f(t );
101
+ % calculate the energy consumption of mobile execution
102
+ E(t , 1 ) = k * W * (f(t )^2 );
103
+ % calculate the value of optimization goal
104
+ J_m = - B_hat(t ) * k * W * (f(t ))^2 + V * W / f(t );
105
+ else
106
+ % the sub-problem is not fasible because (i) the limited
107
+ % computation capacity or (ii) time cosumed out of deadline or
108
+ % (iii) the energy consumed out of battery energy level
109
+ % If it is not feasible, it just means that we cannot choose
110
+ % 'I_m=1'. It dosen't mean that the task has to be dropped.
111
+ disp(' mobile execution ($\mathcal{P}_{ME}$) is not feasible!' )
112
+ f(t ) = 0 ;
113
+ cost(t , 1 ) = 0 ;
114
+ E(t , 1 ) = 0 ;
115
+ % 'I_m=1' can never be chosen if mobile execution goal is inf
116
+ J_m = inf ;
117
+ % Attention! We do not check whether the energy cunsumed is larger than
118
+ % battery energy level because the problem $\mathcal{J}_{CO}$ does
119
+ % not have constraint (8).
120
+ end
121
+
122
+ %% step 2.2: solve the optimization problem $\mathcal{P}_{SE}$ (p(t) > 0)
123
+ E_tmp = sigma * L * log(2 ) / (omega * h );
124
+ p_L_taud = (power(2 , L / (omega * tau_d )) - 1 ) * sigma / h ;
125
+ % calculate p_L
126
+ if E_tmp >= E_min
127
+ p_L = p_L_taud ;
128
+ else
129
+ % calculate p_E_min (use inline function and fsolve)
130
+ y = @(x ) x * L - omega * log2(1 + h * x / sigma ) * E_min ;
131
+ % accroding to the function figure, p_L_taud is a positive
132
+ % number around 0.2
133
+ p_E_min = fsolve(y , 0.2 , opt );
134
+ p_L = max(p_L_taud , p_E_min );
135
+ end
136
+ % calculate p_U
137
+ if E_tmp >= E_max
138
+ p_U = 0 ;
139
+ else
140
+ % caculate p_E_max (use inline function and fsolve)
141
+ y = @(x ) x * L - omega * log2(1 + h * x / sigma ) * E_max ;
142
+ % accroding to the function figure, p_E_max is a large positive
143
+ % number around 20
144
+ p_E_max = fsolve(y , 100 , opt );
145
+ p_U = min(p_tx_max , p_E_max );
146
+ end
147
+ if p_L <= p_U
148
+ % the sub-problem is feasible
149
+ disp(' MEC server execution ($\mathcal{P}_{SE}$) is feasible!' )
150
+ % calculate p_0
151
+ virtual_battery = B_hat(t );
152
+ y = @(x ) virtual_battery * log2(1 + h * x / sigma ) + ...
153
+ h * (V - virtual_battery * x ) / log(2 ) / (sigma + h * x );
154
+ p_0 = fsolve(y , 0.5 , opt );
155
+
156
+ if (p_U < p_0 && B_hat(t ) < 0 ) || B_hat(t ) >= 0
157
+ p(t ) = p_U ;
158
+ elseif p_0 < p_L && B_hat(t ) < 0
159
+ p(t ) = p_L ;
160
+ elseif p_0 >= p_L && p_0 <= p_U && B_hat(t ) < 0
161
+ p(t ) = p_0 ;
162
+ end
163
+ % check whether p(t) is zero
164
+ if p(t ) == 0
165
+ disp(' Something wrong! p is 0!' )
166
+ end
167
+
168
+ % calculate the delay of MEC server execution
169
+ cost(t , 2 ) = L / (omega * log2(1 + h * p(t )/sigma ));
170
+ % calculate the energy consumption of MEC server execution
171
+ E(t , 2 ) = p(t ) * cost(t , 2 );
172
+ % calculate the value of optimization goal
173
+ J_s = (-B_hat(t ) * p(t ) + V ) * cost(t , 2 );
174
+ else
175
+ % the sub-problem is not feasible because (i) the limited transmit
176
+ % power or (ii) time cosumed out of deadline or (iii) the energy
177
+ % consumed out of battery energy level
178
+ % If it is not feasible, it just means that we cannot choose
179
+ % 'I_s=1'. It dosen't mean that the task has to be dropped.
180
+ disp(' MEC server execution ($\mathcal{P}_{SE}$) is not feasible!' )
181
+ p(t ) = 0 ;
182
+ cost(t , 2 ) = 0 ;
183
+ E(t , 2 ) = 0 ;
184
+ % 'I_s=1' can never be chosen if MEC server execution goal is inf
185
+ J_s = inf ;
186
+ % Similarly, we do not check whether the energy cunsumed is larger than
187
+ % battery energy level because the problem $\mathcal{J}_{CO}$ does
188
+ % not have constraint (8).
189
+ end
190
+
191
+ %% step 3: choose the best execution mode
192
+ J_d = V * phi ;
193
+ disp([' J_m:' , num2str(J_m )])
194
+ disp([' J_s:' , num2str(J_s )])
195
+ [~ , mode ] = min([J_m , J_s , J_d ]);
196
+ chosen_mode(t ) = mode ;
197
+ end
198
+
199
+ %% step 4: according to the chosen execution mode, calculate the real dealy and energy consumption
200
+ if chosen_mode(t ) == 1
201
+ % mobile execution is chosen
202
+ cost(t , 3 ) = cost(t , 1 );
203
+ E(t , 3 ) = E(t , 1 );
204
+ elseif chosen_mode(t ) == 2
205
+ % MEC server execution is chosen
206
+ cost(t , 3 ) = cost(t , 2 );
207
+ E(t , 3 ) = E(t , 2 );
208
+ elseif chosen_mode(t ) == 3
209
+ % task is dropped, the delay is the task dropping penalty and the
210
+ % energy consumption is zero
211
+ cost(t , 3 ) = phi ;
212
+ E(t , 3 ) = 0 ;
213
+ else
214
+ % no task is requested, the delay and the energy consumption are
215
+ % both zero
216
+ cost(t , 3 ) = 0 ;
217
+ E(t , 3 ) = 0 ;
218
+ end
219
+
220
+ %% step 5: update the battery energy level and go to next time slot
221
+ B(t + 1 ) = B(t ) - E(t , 3 ) + e(t );
222
+ t = t + 1 ;
223
+ end
224
+
225
+ %% step 6: evaluate the simulation results
226
+ % 1. the battery energy level vs. time slot
227
+ subplot(2 , 1 , 1 );
228
+ plot(1 : T , B(1 : T ));
229
+ hold on
230
+ plot(1 : T , repmat(theta + E_H_max , [T , 1 ]), ' -' )
231
+ title(' Envolution of battery energy level' )
232
+ xlabel(' time slot' )
233
+ ylabel(' battery energy level $B_t$' , ' Interpreter' ,' latex' )
234
+
235
+ % 2. the average execution cost vs. time slot
236
+ accumulated = 0 ;
237
+ average_cost = zeros(T , 1 );
238
+ request_num = 0 ;
239
+ for t = 1 : T
240
+ accumulated = accumulated + cost(t , 3 );
241
+ if cost(t , 3 ) ~= 0
242
+ % there exists task request
243
+ request_num = request_num + 1 ;
244
+ end
245
+ average_cost(t ) = accumulated / request_num ;
246
+ end
247
+ subplot(2 , 1 , 2 );
248
+ plot(1 : T , average_cost );
249
+ title(' Envolution of average execution cost' )
250
+ xlabel(' time slot' )
251
+ ylabel(' average execution cost $\f rac{1}{T} \sum_{t=0}^{T-1} cost^t$' , ' Interpreter' ,' latex' )
0 commit comments