1function [UN,QN,p2opt]=qrf_noblo_mmi_ld(MAPs,N,rt,alpha)
3% f; % finite capacity queue
4% M, integer, > 0; % number of queues
5% MR, integer, > 0; % number of independent blocking configurations
6% MM {m in 1:MR, i in 1:M} >=0; % blocking order
7% MM1 {m in 1:MR, i in 1:M}; % blocking order
8% ZZ {m in 1:MR} >=0; % nonzeros in independent blocking configurations
9% ZM, integer, >=0; % max of ZZ
10% BB {m in 1:MR, i in 1:M} >=0; % blocking state
11% K {i in 1:M}, integer, > 0; % number of phases
for each queue
12% F {i in 1:M}, integer, > 0; % capacity
13% N, integer, >0; % population
14% mu {i in 1:M, k in 1:K(i), h in 1:K(i)} >=0; % completion transition rates
15% v {i in 1:M, k in 1:K(i), h in 1:K(i)} >=0; % background transition rates
16% r {i in 1:M, j in 1:M} >=0; % routing probabilities
19%var p2 {j = 1:M, nj = 1+(0:N), k = 1:K(j), i = 1:M, ni = 1+(0:N), h = 1:K(i), m = 1:MR} >= 0;
20%var e {i = 1:M, k = 1:K(i)} >=0;
32 K(i)=size(MAPs{i}{1},1);
43 for h=1:size(MAPs{i}{1},1)
44 for k=1:size(MAPs{i}{1},1)
45 mu(i,h,k)=MAPs{i}{2}(h,k);
51 for h=1:size(MAPs{i}{1},1)
52 for k=1:size(MAPs{i}{1},1)
56 v(i,k,h)=MAPs{i}{1}(h,k);
62q = zeros(M,M,max(K),max(K));
70 q(i,j,k,h,1+ni,1+nj) = rt(i,j)*mu(i,k,h)*alpha(i,ni);
72 q(i,j,k,h,1+ni,1+nj) = v(i,k,h)+rt(i,i)*mu(i,k,h)*alpha(i,ni);
81px = ones(M*max(K),1); px = px/sum(px(:));
82x = [ones(M*(N+1)*max(K)*M*(N+1)*max(K)*MR,1); px];
85options = optimset(
'fmincon');
86options.Display =
'iter';
87%options.LargeScale =
'off';
89options.MaxFunEvals = 1e3;
90options.MaxSQPIter = 50;
92options.Algorithm =
'sqp';
93%options.OutputFcn = @outfun;
95[xopt, fopt] = fmincon(@(x) mmi(x),x,[],[],[],[],x*0,x*0+1,@(x) sub_qrfcon(x,q,M,MR,MM,MM1,ZZ,ZM,BB,F,N),options);
96[p2opt,~] = sub_qrfvar(xopt);
104 UN(ti) = UN(ti) + p2opt(ti,ni,ki,ti,ni,ki, m);
105 QN(ti) = QN(ti) + ni*p2opt(ti,ni,ki,ti,ni,ki, m);
111% add LB>=0 to both p2 and e
112% LINEAR PROGRAMMING: UTILIZATION UPPER BOUND AT QUEUE 1
113%minimize U1min: sum {m in 1..MR} sum {k in 1..K[1]} sum {n1 in 1..F[1]} p2[1,n1,k,1,n1,k,m];
116 function fobj = mmi(x)
118 %minimize MI: sum {m in 1..MR} sum {i in 1..M, j in 1..M, ki in 1..K[i], kj in 1..K[j]: i<>j} sum {nj in 0..F[j]} sum {ni in 0..F[j]} p2[i,ni,ki,j,nj,kj,m]*(log(1e-6+p2[i,ni,ki,j,nj,kj,m])-log(1e-6+p2[i,ni,ki,i,ni,ki,m])-log(1e-6+p2[j,nj,kj,j,nj,kj,m]));
119 [p2,~] = sub_qrfvar(x);
129 fobj = fobj + p2(i,ni,ki,j,nj,kj,m)*(log(1e-6+p2(i,ni,ki,j,nj,kj,m))-log(1e-6+p2(i,ni,ki,i,ni,ki,m))-log(1e-6+p2(j,nj,kj,j,nj,kj,m)));
140 function fobj = mem(x)
142 %maximize H: -sum {m in 1..MR} sum {i in 1..M} sum {k in 1..K[i]} sum {ni in 1..F[i]} p2[i,ni,k,i,ni,k,m]*log(1e-6+p2[i,ni,k,i,ni,k,m]);
143 [p2,~] = sub_qrfvar(x);
149 fobj = fobj - p2(i,ni,k,i,ni,k,m)*log(1e-6 + p2(i,ni,k,i,ni,k,m));
156 function fobj = umin(x, i)
157 [p2,~] = sub_qrfvar(x);
158 fobj = sum(p2(i,1+(1:F(i)),1:K(i),i,1+(1:F(i)),1:K(i), 1:MR),
'all');
161 function [p2,e] = sub_qrfvar(x)
163 p2 = zeros(M,N+1,max(K),M,N+1,max(K),MR);
171 p2(j,nj,k,i,ni,h,m) = x(ctr);
189 function [c,ceq] = sub_qrfcon(x,q,M,MR,MM,MM1,ZZ,ZM,BB,F,N)
194 [p2,e] = sub_qrfvar(x);
197 % subject to ONE {j in 1..M}: sum {nj in 0..N, k in 1..K[j], m in 1..MR} p2[j,nj,k,j,nj,k,m]=1;
201 for nj = 1+(0:N),
for k = 1:K(j), for m = 1:MR
202 ceq(end) = ceq(end) + p2(j,nj,k,j,nj,k,m);
205 ceq(end) = ceq(end) -1;
208 % subject to ZERO1 {j in 1..M, k in 1..K[j], nj in 0..N, i in 1..M, h in 1..K[i], ni in 0..N, m in 1..MR: i==j and nj==ni and h<>k}: p2[j,nj,k,i,ni,h,m]=0;
209 for j = 1:M,
for k =1:K(j),
for nj = 1+(0:N),
for i = 1:M, for h = 1:K(i), for ni = 1+(0:N), for m = 1:MR
210 if i==j && (nj-1)==(ni-1) && h~=k % rescaled back nj and ni
211 ceq(end+1) = p2(j,nj,k,i,ni,h,m); %=0
213 end, end, end, end, end, end, end
215 % subject to ZERO2 {j in 1..M, k in 1..K[j], nj in 0..N, i in 1..M, h in 1..K[i], ni in 0..N, m in 1..MR: i==j and nj<>ni}: p2[j,nj,k,i,ni,h,m]=0;
216 for j = 1:M,
for k =1:K(j),
for nj = 1+(0:N),
for i = 1:M, for h = 1:K(i), for ni = 1+(0:N), for m = 1:MR
217 if i==j && (nj-1)~=(ni-1) % rescaled back nj and ni
218 ceq(end+1) = p2(j,nj,k,i,ni,h,m); %=0
220 end, end, end, end, end, end, end
222 % subject to ZERO3 {j in 1..M, k in 1..K[j], nj in 0..N, i in 1..M, h in 1..K[i], ni in 0..N, m in 1..MR: i<>j and nj+ni>N}: p2[j,nj,k,i,ni,h,m]=0;
223 for j = 1:M,
for k =1:K(j),
for nj = 1+(0:N),
for i = 1:M, for h = 1:K(i), for ni = 1+(0:N), for m = 1:MR
224 if i~=j && (nj-1)+(ni-1)>N % rescaled back nj and ni
225 ceq(end+1) = p2(j,nj,k,i,ni,h,m);
227 end, end, end, end, end, end, end
229 % subject to ZERO5 {j in 1..M, k in 1..K[j], i in 1..M, h in 1..K[i], ni in 0..F[i], m in 2..MR: BB[m,j]==1}: p2[j,0,k,i,ni,h,m]=0;
230 for j = 1:M,
for k =1:K(j),
for i = 1:M,
for h = 1:K(i),
for ni = 1+(0:F(i)),
for m = 2:MR
232 ceq(end+1) = p2(j,1+0,k,i,ni,h,m);
234 end, end, end, end, end, end
236 % subject to ZERO6 {j in 1..M, k in 1..K[j], nj in F[j]+1..N, i in 1..M, h in 1..K[i], ni in 0..N, m in 1..MR}: p2[j,nj,k,i,ni,h,m]=0;
237 for j = 1:M,
for k =1:K(j),
for nj = 1+((F(j)+1):N), for i = 1:M, for h = 1:K(i), for ni = 1+(0:N), for m = 1:MR
238 ceq(end+1) = p2(j,nj,k,i,ni,h,m);
239 end, end, end, end, end, end, end
241 % subject to ZERO7 {j in 1..M, k in 1..K[j], nj in 1..F[j], i in 1..M, h in 1..K[i], ni in 0..N, m in 2..MR: BB[m,j]==1 and i<>j and i<>f and ni+nj+F[f]>N}: p2[j,nj,k,i,ni,h,m]=0;
242 for j = 1:M,
for k =1:K(j),
for nj = 1+(1:F(j)),
for i = 1:M, for h = 1:K(i), for ni = 1+(0:N), for m = 2:MR
243 if BB(m,j)==1 && i~=j && i~=f && (ni-1)+(nj-1)+F(f)>N % rescaled back ni and nj
244 ceq(end+1) = p2(j,nj,k,i,ni,h,m);
246 end, end, end, end, end, end, end
248 % subject to SIMMETRY {j in 1..M, nj in 0..N, k in 1..K[j], i in 1..M, ni in 0..N, h in 1..K[i], m in 1..MR}: p2[i,ni,h,j,nj,k,m] = p2[j,nj,k,i,ni,h,m];
249 for j = 1:M,
for nj = 1+(0:N),
for k =1:K(j), for i = 1:M, for ni = 1+(0:N), for h = 1:K(i), for m = 1:MR
250 ceq(end+1) = p2(i,ni,h,j,nj,k,m) - p2(j,nj,k,i,ni,h,m);
251 end, end, end, end, end, end, end
253 % subject to MARGINALS {j in 1..M, k in 1..K[j], nj in 0..N, i in 1..M, m in 1..MR: i<>j}: p2[j,nj,k,j,nj,k,m]= sum {ni in 0..N-nj} sum {h in 1..K[i]} p2[j,nj,k,i,ni,h,m];
254 for j = 1:M,
for k =1:K(j),
for nj = 1+(0:N),
for i = 1:M, for m = 1:MR
257 ceq(end+1) = p2(j,nj,k,j,nj,k,m);
259 for ni = 1+(0:(N-nj)) %% added +1
261 ceq(end) = ceq(end) - p2(j,nj,k,i,ni,h,m);
265 end, end, end, end, end
267 %subject to UEFF {j in 1..M, i in 1..M, ki in 1..K[i]}: e[i,ki] = sum {nj in 0..N, kj in 1..K[j], m in 1..MR, ni in 1..N: BB[m,i]==0} p2[j,nj,kj,i,ni,ki,m];
268 for j = 1:M,
for i = 1:M,
for ki = 1:K(i)
270 ceq(end+1) = e(i,ki);
272 for nj = 1+(0:N),
for kj = 1:K(j), for m = 1:MR, for ni = 1+(1:N)
274 ceq(end) = ceq(end) - p2(j,nj,kj,i,ni,ki,m);
279 %subject to THM1 {i in 1..M, k in 1..K[i]}: sum {j in 1..M, h in 1..K[i]} q[i,j,k,h]*e[i,k] =sum {j in 1..M, h in 1..K[i]} q[i,j,h,k]*e[i,h];
280 for i = 1:M,
for k =1:K(i)
283 for j = 1:M,
for h = 1:K(i)
284 ceq(end) = ceq(end) + q(i,j,k,h)*e(i,k);
287 for j = 1:M,
for h = 1:K(i)
288 ceq(end) = ceq(end) - q(i,j,h,k)*e(i,h);
292 %subject to THM2 {j in 1..M, k in 1..K[j], nj in 0..F[j], m in 1..MR}: sum {i in 1..M, ni in 1..F[i], ki in 1..K[i]} ni*p2[j,nj,k,i,ni,ki,m]= N*p2[j,nj,k,j,nj,k,m];
293 for j = 1:M,
for k =1:K(j),
for nj = 1+(0:F(j)),
for m = 1:MR
296 for i = 1:M,
for ni = 1+(1:F(i)),
for ki = 1:K(i)
297 ceq(end) = ceq(end) + (ni-1)*p2(j,nj,k,i,ni,ki,m); % recaled back ni
300 ceq(end) = ceq(end) - N*p2(j,nj,k,j,nj,k,m);
303 %subject to COR1 : sum {m in 1..MR, i in 1..M, j in 1..M, nj in 1..F[j], ni in 1..F[i], ki in 1..K[i], kj in 1..K[j]} ni*nj*p2[j,nj,kj,i,ni,ki,m]= N^2;
306 for m = 1:MR,
for i = 1:M,
for j = 1:M,
for nj = 1+(1:F(j)),
for ni = 1+(1:F(i)), for ki = 1:K(i), for kj = 1:K(j)
307 ceq(end) = ceq(end) + (ni-1)*(nj-1)*p2(j,nj,kj,i,ni,ki,m); % rescaled back ni and nj
308 end, end, end, end, end, end, end
310 ceq(end) = ceq(end) - N^2;
312 %subject to THM30 {i in 1..M, u in 1..K[i]: i<>f}: sum {j in 1..M, nj in 1..F[j], k in 1..K[j], h in 1..K[j], m in 1..MR: j<>i and j<>f and BB[m,j]==0}
313 % q[j,i,k,h]*p2[j,nj,k, i,0,u, m] + sum {j in 1..M, nj in 1..F[j], k in 1..K[j], h in 1..K[j], m in 1..MR: j<>i and j==f and MM[m,1]<>i} q[j,i,k,h]*p2[j,nj,k, i,0,u, m]
314 % = sum {j in 1..M, nj in 0..F[j], k in 1..K[i], h in 1..K[j], m in 1..MR: j<>i and j<>f and BB[m,i]==0} q[i,j,k,u]*p2[j,nj,h, i,0+1,k, m]
315 % + sum {j in 1..M, nj in 0..F[j], k in 1..K[i], h in 1..K[j], m in 1..MR: j<>i and j==f and BB[m,i]==0 and nj<F[j]} q[i,j,k,u]*p2[j,nj,h, i,1,k, m]
316 % + sum {j in 1..M, nj in 0..F[j], y in 1..K[j], m in 1..MR: j<>i and j==f and BB[m,i]==1 and nj==F[j] and MM[m,1]==i} sum {p in 1..K[f], w in 1..M: w<>f and w<>i} q[f,w,y,p]*p2[f,nj,y, i,1,u, m] ;
318 f = -1; % never equal to FC station
319 for i = 1:M,
for u = 1:K(i)
323 for j = 1:M,
for nj = 1+(1:F(j)),
for k = 1:K(j), for h = 1:K(j), for m = 1:MR
324 if j~=i && j~=f && BB(m,j)==0
325 ceq(end)= ceq(end) + q(j,i,k,h)*p2(j,nj,k, i,1+0,u, m);
327 end, end, end, end, end
329 for j = 1:M,
for nj = 1+(1:F(j)),
for k = 1:K(j), for h = 1:K(j), for m = 1:MR
330 if j~=i && j==f && MM(m,1)~=i
331 ceq(end)= ceq(end) + q(j,i,k,h)*p2(j,nj,k, i,1+0,u, m);
333 end, end, end, end, end
335 for j = 1:M,
for nj = 1+(0:F(j)),
for k = 1:K(i), for h = 1:K(j), for m = 1:MR
336 if j~=i && j~=f && BB(m,i)==0
337 ceq(end)= ceq(end) - q(i,j,k,u)*p2(j,nj,h, i,1+0,k, m);
339 end, end, end, end, end
341 for j = 1:M,
for nj = 1+(0:F(j)),
for k = 1:K(i), for h = 1:K(j), for m = 1:MR
342 if j~=i && j==f && BB(m,i)==0 && (nj-1)<F(j) % rescaled back nj
343 ceq(end)= ceq(end) - q(i,j,k,u)*p2(j,nj,h, i,1+1,k, m);
345 end, end, end, end, end
347 for j = 1:M,
for nj = 1+(0:F(j)),
for y = 1:K(j), for m = 1:MR
348 if j~=i && j==f && BB(m,i)==1 && (nj-1)==F(j) && MM(m,1)==i % rescaled back nj
349 for p = 1:K(f), for w = 1:M
351 ceq(end)= ceq(end) - q(f,w,y,p)*p2(f,nj,y, i,1+1,u, m);
359 % subject to THM3 {i in 1..M, ni in 0..(F[i]-1): i<>f}:
360 % sum {j in 1..M, nj in 1..F[j], k in 1..K[j], h in 1..K[j], u in 1..K[i], m in 1..MR: j<>i and j<>f and BB[m,j]==0} q[j,i,k,h]*p2[j,nj,k, i,ni,u, m]
361 % + sum {j in 1..M, nj in 1..F[j], k in 1..K[j], h in 1..K[j], u in 1..K[i], m in 1..MR: j<>i and j==f and MM[m,1]<>i} q[j,i,k,h]*p2[j,nj,k, i,ni,u, m]
362 % = sum {j in 1..M, nj in 0..F[j], k in 1..K[i], u in 1..K[j], m in 1..MR: j<>i and j<>f and BB[m,i]==0} sum {h in 1..K[i]} q[i,j,k,h]*p2[j,nj,u, i,ni+1,k, m]
363 % + sum {j in 1..M, nj in 0..F[j], k in 1..K[i], u in 1..K[j], m in 1..MR: j<>i and j==f and BB[m,i]==0 and nj<F[j]} sum {h in 1..K[i]} q[i,j,k,h]*p2[j,nj,u, i,ni+1,k, m]
364 % + sum {j in 1..M, nj in 0..F[j], k in 1..K[i], u in 1..K[j], m in 1..MR: j<>i and j==f and BB[m,i]==1 and nj==F[j] and MM[m,1]==i} sum {p in 1..K[f], w in 1..M: w<>f and w<>i} q[f,w,u,p]*p2[f,nj,u, i,ni+1,k, m] ;
365 for i = 1:M,
for ni = 1+(0:(F(i)-1))
369 for j = 1:M,
for nj = 1+(1:F(j)),
for k = 1:K(j), for h = 1:K(j), for u = 1:K(i), for m = 1:MR
370 if j~=i && j~=f && BB(m,j)==0
371 ceq(end)= ceq(end) + q(j,i,k,h)*p2(j,nj,k, i,ni,u, m);
373 end, end, end, end, end, end
374 for j = 1:M,
for nj = 1+(1:F(j)),
for k = 1:K(j), for h = 1:K(j), for u = 1:K(i), for m = 1:MR
375 if j~=i && j==f && MM(m,1)~=i
376 ceq(end)= ceq(end) + q(j,i,k,h)*p2(j,nj,k, i,ni,u, m);
378 end, end, end, end, end, end
380 for j = 1:M,
for nj = 1+(0:F(j)),
for k = 1:K(i), for u = 1:K(j), for m = 1:MR
381 if j~=i && j~=f && BB(m,i)==0
383 ceq(end)= ceq(end) - q(i,j,k,h)*p2(j,nj,u, i,ni+1,k, m);
386 end, end, end, end, end
387 for j = 1:M,
for nj = 1+(0:F(j)),
for k = 1:K(i), for u = 1:K(j), for m = 1:MR
388 if j~=i && j==f && BB(m,i)==0 && (nj-1)<F(j) % rescaled back nj
390 ceq(end)= ceq(end) - q(i,j,k,h)*p2(j,nj,u, i,ni+1,k, m);
393 end, end, end, end, end
394 for j = 1:M,
for nj = 1+(0:F(j)),
for k = 1:K(i), for u = 1:K(j), for m = 1:MR
395 if j~=i && j==f && BB(m,i)==1 && (nj-1)==F(j) && MM(m,1)==i % rescaled back nj
396 for p = 1:K(f), for w = 1:M
398 ceq(end)= ceq(end) - q(f,w,u,p)*p2(f,nj,u, i,ni+1,k, m);
402 end, end, end, end, end
406 % subject to THM3f {i in 1..M, ni in 0..(F[i]-1): i==f}: sum {j in 1..M, nj in 1..F[j], k in 1..K[j], h in 1..K[j], u in 1..K[i], m in 1..MR: j<>i and j<>f and BB[m,j]==0 and ni < F[i]} q[j,i,k,h]*p2[j,nj,k, i,ni,u, m]
407 % + sum {j in 1..M, nj in 1..F[j], k in 1..K[j], h in 1..K[j], u in 1..K[i], m in 1..MR: j<>i and j==f and ni==F[i] and MM[m,1]==j} sum {w in 1..M: w<>f} q[j,i,k,h]*p2[j,nj,k, i,ni,u, m]
408 % = sum {j in 1..M, nj in 0..F[j], k in 1..K[i], h in 1..K[i], u in 1..K[j]: j<>i and ni < F[i]} q[i,j,k,h]*p2[j,nj,u, i,ni+1,k, 1];
409 for i = 1:M,
for ni = 1+(0:(F(i)-1))
413 for j = 1:M,
for nj = 1+(1:F(j)),
for k = 1:K(j), for h = 1:K(j), for u = 1:K(i), for m = 1:MR
414 if j~=i && j~=f && BB(m,j)==0 && (ni-1) < F(i) % rescaled back ni
415 ceq(end) = ceq(end) + q(j,i,k,h)*p2(j,nj,k, i,ni,u, m);
417 end, end, end, end, end, end
418 for j = 1:M,
for nj = 1+(1:F(j)),
for k = 1:K(j), for h = 1:K(j), for u = 1:K(i), for m = 1:MR
419 if j~=i && j==f && (ni-1)==F(i) && MM(m,1)==j % rescaled back ni
422 ceq(end) = ceq(end) + q(j,i,k,h)*p2(j,nj,k, i,ni,u, m);
426 end, end, end, end, end, end
428 for j = 1:M,
for nj = 1+(0:F(j)),
for k = 1:K(i), for h = 1:K(i), for u = 1:K(j)
429 if j~=i && (ni-1) < F(i) % rescaled back ni
430 ceq(end) = ceq(end) - q(i,j,k,h)*p2(j,nj,u, i,ni+1,k, 1);
432 end, end, end, end, end
436 % subject to THM4 {j in 1..M, k in 1..K[j], i in 1..M, m in 1..MR}: sum{t in 1..M} sum {h in 1..K[t]} sum {nj in 0..N} sum {nt in 0..N} nt*p2[j,nj,k,t,nt,h,m]
437 % >= N*sum {h in 1..K[i]} sum {nj in 0..N} sum {ni in 1..N} (p2[j,nj,k,i,ni,h,m]);
438 for j = 1:M,
for k = 1:K(j),
for i = 1:M,
for m = 1:MR
439 c(end+1) = 0; % <= inequality
440 % LHS with sign swapped since >= in GLPK
441 for t = 1:M,
for h = 1:K(t),
for nj = 1+(0:N),
for nt = 1+(0:N)
442 c(end) = c(end) - (nt-1)*p2(j,nj,k,t,nt,h,m); % rescaled back nt
444 % RHS with sign swapped since >= in GLPK
445 for h = 1:K(i),
for nj = 1+(0:N),
for ni = 1+(1:N)
446 c(end) = c(end) + N*(p2(j,nj,k,i,ni,h,m));