-
Notifications
You must be signed in to change notification settings - Fork 3
/
Copy pathOptimiseDistributed2Beta.m
188 lines (157 loc) · 6.8 KB
/
OptimiseDistributed2Beta.m
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
185
186
187
188
%% Test effect of Beta *after DeleteDistantStations
% Beta is a parameter in the distributed algorithm
%clearvars -except Arrival capacity e eta iifl NbhDistance NumStation
addpath(genpath('/Applications/MATLAB_R2014b.app/toolbox/yalmip/'));
% add yalmip solver to matlab search path
tic
%% Initialise optimisation variables/constraints/objective
global beta k_stack
k_stack = Inf(72,1);
%ProblemSize = NumStation;
ProblemSize = 100;
distance = NbhDistance(1:ProblemSize,1:ProblemSize);
NbhVec = e(1:ProblemSize,1:ProblemSize);
Deviation_Indicator = 0;
solution = Inf(ProblemSize,ProblemSize,72); % decision variable
for Tslice = 1:71
%% Fill level dynamics
ticTslice = tic;
if Tslice == 1
fl = iifl(1:ProblemSize);
else
fl = fl + eta(1:ProblemSize,Tslice-1) ...
+ sum(solution(:,:,Tslice-1),1)' ...
- sum(solution(:,:,Tslice-1),2);
% carry the previous fill level
% include the net change eta
% consider the extra #bikes deviated to here at time t-1
% consider #bikes left at time t-1
end
emptylevel = capacity(1:ProblemSize)-fl;
fl(fl<0) = 0; % under-empty stations will be set as empty
% some customers would not be able to depart
% It should not be over-full after optimising
if min(emptylevel) < 0 % If so, error message
fprintf('Station %d with empty level %d at Tslice %d\n', ...
find(emptylevel == min(emptylevel),1), min(emptylevel), Tslice);
end
%% Distributed Constrained Optimisation
skipflag = false(1);
converged = false(1);
%beta = 20;
c = Inf(ProblemSize,1);
k_max = 50;
k_si = k_max;
k = 1;
primal_obj_x = Inf(k_max,1);
primal_obj_xhat = Inf(k_max,1);
viola_x = Inf(k_max,1);
viola_xhat = Inf(k_max,1);
x = Inf(ProblemSize,ProblemSize,k_max); % row vectors
xhat = Inf(ProblemSize,ProblemSize,k_max); % row vectors
xhat(:,:,1) = zeros(ProblemSize,ProblemSize);
for i = 1:ProblemSize
xhat(i,i,1) = Arrival(i,Tslice);
% The optimised case neglecting capacity-fl inequality constraints
end
lamb = Inf(ProblemSize,ProblemSize,k_max); % row vectors
lamb(:,:,1) = zeros(ProblemSize,ProblemSize); % lamb(1) = 0
l = Inf(ProblemSize,ProblemSize,k_max); % row vectors
lamb_convergence = Inf(k_max,1); % lamb convergence rate at each iteration
% Skip computation if no Arrivals at all
if sum(Arrival(1:ProblemSize,Tslice)) == 0
skipflag = true(1);
solution(:,:,Tslice) = zeros(ProblemSize,ProblemSize);
end
%% Repeat until convergence
sum_tocI = zeros(ProblemSize,1);
if ~skipflag % if there are some Arrivals
while k<k_max && ~converged
c(k) = beta/k;
% Implement two different sequences for xhat/xtilta
if k < k_si
cc = c(k)/sum(c(1:k));
else
cc = c(k)/sum(c(k_si:k));
end
diagonal = min(Arrival(1:ProblemSize,Tslice),emptylevel);
fulllist = zeros(ProblemSize,1);
countfull = 0;
for checki = 1:ProblemSize
if emptylevel(checki) - diagonal(checki) == 0
countfull = countfull+1;
fulllist(countfull) = checki;
end
end
fulllist = fulllist(1:countfull);
for i = 1:ProblemSize
ticI = tic;
l(i,:,k) = mean(lamb(:,:,k),1); % by considering a(i,j) = 1/m
opt = sdpvar(ProblemSize,1); % row vector illustrated as a column
constraints = [opt >= 0, sum(opt) == Arrival(i,Tslice), opt(i) == diagonal(i)];
for counti = 1:countfull % if a station is full, do not receive any
if fulllist(counti) ~= i
constraints = [constraints, opt(fulllist(counti)) == 0];
end
end
%local constraints
objective = distance(i,:) * opt + l(i,:,k) * (NbhVec(:,i) .* opt - emptylevel./ProblemSize);
options = sdpsettings('verbose',0,'solver','linprog');
sol = optimize(constraints,objective,options);
if sol.problem == 0 % no problem
x(i,:,k+1) = value(opt); % row vector
primal_obj_x(k+1) = value(objective);
else
display('Something went wrong');
sol.info
yalmiperror(sol.problem)
end
lamb_next = l(i,:,k) + c(k) * (x(i,:,k+1) - emptylevel'./ProblemSize);
lamb_next(lamb_next <0) = 0; %projection
lamb(i,:,k+1) = lamb_next;
xhat(i,:,k+1) = xhat(i,:,k) + cc * (x(i,:,k+1)-xhat(i,:,k));
tocI = toc(ticI);
%fprintf('Agent %d with elapsed time', i);toc(ticI);
sum_tocI(i) = sum_tocI(i) + tocI;
end
primal_obj_x(k+1) = sum(x(:,:,k+1) * distance(:,i));
primal_obj_xhat(k+1) = sum(xhat(:,:,k+1) * distance(:,i));
v = sum(x(:,:,k+1),1) - emptylevel'; v(v <0) = 0; viola_x(k+1) = sum(v);
v = sum(xhat(:,:,k+1),1) - emptylevel'; v(v <0) = 0; viola_xhat(k+1) = sum(v);
%%% Check Convergence
if k > 1
temp = 0;
for checki = 1:ProblemSize
lamb_change = lamb(checki,:,k+1)-lamb(checki,:,k);
temp = temp + norm(lamb_change) ./ norm(lamb(checki,:,k));
end
lamb_convergence(k) = temp ./ ProblemSize;
end
% if almost converges AND still in phase 1
if lamb_convergence(k) < 0.005 && k < k_si
k_si = k+1;
end
converged = true(1);
for checki = 1:ProblemSize
if sum(round(xhat(:,checki,k+1))) > emptylevel(checki) ...
|| sum(round(xhat(checki,:,k+1))) ~= Arrival(checki,Tslice)
% if any column violates filllevel constraints
% or if any row (after round off) gets incorrect #Arrival
converged = false(1);
break;
end
end
k = k + 1;
end
k_end = k;
k_stack(Tslice) = k_end-1;
% Store final solution for time at Tslice
solution(:,:,Tslice) = round(xhat(:,:,k_end));
end
Deviation_Indicator = Deviation_Indicator + sum(sum(solution(:,:,Tslice))) - trace(solution(:,:,Tslice));
if mod(Tslice,10) == 0
fprintf('Time slice computed is %d\n', Tslice);
end
tocTslice = toc(ticTslice);
end
toc