LINE Solver
MATLAB API documentation
Loading...
Searching...
No Matches
qrf_noblo_mmi_linear.m
1function [UN,QN,p2opt]=qrf_noblo_mmi_linear(MAPs,N,rt,alpha)
2%%% PARAMETERS %%%
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
17
18%%% VARIABLES %%%
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;
21
22M = length(MAPs);
23MR = 1;
24MM = 0;
25F = repmat(N,M,1);
26MM1 = ones(1,3);
27ZZ = zeros(1,M);
28ZM = 0;
29BB = zeros(MR,M);
30
31for i=1:M
32 K(i)=size(MAPs{i}{1},1);
33end
34
35r = zeros(M);
36for i=1:M
37 for j=1:M
38 r(i,j)=rt(i,j);
39 end
40end
41
42for i=1:M
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);
46 end
47 end
48end
49
50for i=1:M
51 for h=1:size(MAPs{i}{1},1)
52 for k=1:size(MAPs{i}{1},1)
53 if h==k
54 v(i,k,h)=0;
55 else
56 v(i,k,h)=MAPs{i}{1}(h,k);
57 end
58 end
59 end
60end
61
62q = zeros(M,M,max(K),max(K));
63for i = 1:M
64 for j = 1:M
65 for k = 1:K(i)
66 for h = 1:K(i)
67 for ni=1:N
68 for nj=1:N
69 if j ~= i
70 q(i,j,k,h,1+ni,1+nj) = rt(i,j)*mu(i,k,h)*alpha(i,ni);
71 else
72 q(i,j,k,h,1+ni,1+nj) = v(i,k,h)+rt(i,i)*mu(i,k,h)*alpha(i,ni);
73 end
74 end
75 end
76 end
77 end
78 end
79end
80
81x0 = [ones(M*(N+1)*max(K)*M*(N+1)*max(K)*MR,1); ones(M*max(K),1)];
82x0 = x0 / sum(x0(:));
83
84options = optimset('fmincon');
85options.Display = 'iter';
86options.LargeScale = 'off';
87options.MaxIter = 100;
88options.MaxFunEvals = 1e5;
89options.MaxSQPIter = 500;
90options.TolCon = 1e-8;
91options.Algorithm = 'sqp';
92%options.Algorithm = 'active-set';
93%options.OutputFcn = @outfun;
94
95[A,b,Aeq,beq] = sub_qrfcon(q,M,MR,MM,BB,F,N);
96
97%[xopt, fopt] = fmincon(@(x) mmi(x),x0,A,b,Aeq,beq,x0*0,x0*0+1,[],options);
98[xopt, fopt] = fmincon(@(x) mem(x),x0,A,b,Aeq,beq,x0*0,x0*0+1,[],options);
99[p2opt,~] = sub_qrfvar(xopt);
100
101for ti=1:M
102 UN(ti) = 0;
103 QN(ti) = 0;
104 for m=1:MR
105 for ni=1+(1:F(ti))
106 for ki=1:K(ti)
107 UN(ti) = UN(ti) + p2opt(ti,ni,ki,ti,ni,ki, m);
108 QN(ti) = QN(ti) + ni*p2opt(ti,ni,ki,ti,ni,ki, m);
109 end
110 end
111 end
112end
113
114% add LB>=0 to both p2 and e
115% LINEAR PROGRAMMING: UTILIZATION UPPER BOUND AT QUEUE 1
116%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];
117
118
119 function fobj = mmi(x)
120 % MMI
121 %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]));
122 [p2,~] = sub_qrfvar(x);
123 fobj = 0;
124 for m = 1:MR
125 for i = 1:M
126 for ki = 1:K(i)
127 for j = 1:M
128 if i~=j
129 for kj = 1:K(j)
130 for ni = 1+(1:F(i))
131 for nj = 1+(1:F(j))
132 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)));
133 end
134 end
135 end
136 end
137 end
138 end
139 end
140 end
141 end
142
143 function fobj = mem(x)
144 % MEM
145 %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]);
146 [p2,~] = sub_qrfvar(x);
147 fobj = 0;
148 for m = 1:MR
149 for i = 1:M
150 for k = 1:K(i)
151 for ni = 1+(1:F(i))
152 fobj = fobj - p2(i,ni,k,i,ni,k,m)*log(1e-6 + p2(i,ni,k,i,ni,k,m));
153 end
154 end
155 end
156 end
157 end
158
159 function fobj = umin(x, i)
160 [p2,~] = sub_qrfvar(x);
161 fobj = sum(p2(i,1+(1:F(i)),1:K(i),i,1+(1:F(i)),1:K(i), 1:MR),'all');
162 end
163
164 function [p2,e] = sub_qrfvar(x)
165 ctr = 1;
166 p2 = zeros(M,N+1,max(K),M,N+1,max(K),MR);
167 for j = 1:M
168 for nj = 1+(0:N)
169 for k = 1:K(j)
170 for i = 1:M
171 for ni = 1+(0:N)
172 for h = 1:K(i)
173 for m = 1:MR
174 p2(j,nj,k,i,ni,h,m) = x(ctr);
175 ctr = ctr + 1;
176 end
177 end
178 end
179 end
180 end
181 end
182 end
183 e = zeros(M,max(K));
184 for i=1:M
185 for k=1:K(i)
186 e(i,k) = x(ctr);
187 ctr = ctr + 1;
188 end
189 end
190 end
191
192 function vec = deltae(j,k)
193 % p2 = zeros(M,N+1,max(K),M,N+1,max(K),MR);
194 % e = zeros(M,max(K));
195 vec = sparse(1,M*(N+1)*max(K)*M*(N+1)*max(K)*MR + M*max(K));
196 idx = M*(N+1)*max(K)*M*(N+1)*max(K)*MR;
197 idx = idx + max(K)*(j-1);
198 idx = idx + k;
199 vec(idx) = 1;
200 end
201
202
203 function vec = deltap2(j,nj,k,i,ni,h,m)
204 % p2 = zeros(M,N+1,max(K),M,N+1,max(K),MR);
205 % e = zeros(M,max(K));
206 vec = sparse(1,M*(N+1)*max(K)*M*(N+1)*max(K)*MR + M*max(K));
207 idx = (N+1)*max(K)*M*(N+1)*max(K)*MR*(j-1);
208 idx = idx + max(K)*M*(N+1)*max(K)*MR*(nj-1);
209 idx = idx + M*(N+1)*max(K)*MR*(k-1);
210 idx = idx + (N+1)*max(K)*MR*(i-1);
211 idx = idx + max(K)*MR*(ni-1);
212 idx = idx + MR*(h-1);
213 idx = idx + m;
214 vec(idx) = 1;
215 end
216
217 function [A,b,Aeq,beq] = sub_qrfcon(q,M,MR,MM,BB,F,N)
218 A=sparse(0,M*(N+1)*max(K)*M*(N+1)*max(K)*MR + M*max(K));
219 Aeq=sparse(0,M*(N+1)*max(K)*M*(N+1)*max(K)*MR + M*max(K));
220 b=sparse(0,1);
221 beq=sparse(0,1);
222
223 %% DEFINITIONS
224 % 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;
225 for j = 1:M
226 % LHS
227 Aeq(end+1,:) = 0;
228 for nj = 1+(0:N), for k = 1:K(j), for m = 1:MR
229 Aeq(end,:) = Aeq(end,:) + deltap2(j,nj,k,j,nj,k,m);
230 end, end, end
231 beq(end+1,1) = -1;
232 % RHS
233 end
234
235 % 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;
236 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
237 if i==j && (nj-1)==(ni-1) && h~=k % rescaled back nj and ni
238 Aeq(end+1,:) = deltap2(j,nj,k,i,ni,h,m); %=0
239 beq(end+1,1) = 0;
240 end
241 end, end, end, end, end, end, end
242
243 % 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;
244 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
245 if i==j && (nj-1)~=(ni-1) % rescaled back nj and ni
246 Aeq(end+1,:) = deltap2(j,nj,k,i,ni,h,m); %=0
247 beq(end+1,1) = 0;
248 end
249 end, end, end, end, end, end, end
250
251 % 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;
252 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
253 if i~=j && (nj-1)+(ni-1)>N % rescaled back nj and ni
254 Aeq(end+1,:) = deltap2(j,nj,k,i,ni,h,m);
255 beq(end+1,1) = 0;
256 end
257 end, end, end, end, end, end, end
258
259 % 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;
260 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
261 if BB(m,j)==1
262 Aeq(end+1,:) = deltap2(j,1+0,k,i,ni,h,m);
263 beq(end+1,1) = 0;
264 end
265 end, end, end, end, end, end
266
267 % 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;
268 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
269 Aeq(end+1,:) = deltap2(j,nj,k,i,ni,h,m);
270 beq(end+1,1) = 0;
271 end, end, end, end, end, end, end
272
273 % 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;
274 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
275 if BB(m,j)==1 && i~=j && i~=f && (ni-1)+(nj-1)+F(f)>N % rescaled back ni and nj
276 Aeq(end+1,:) = deltap2(j,nj,k,i,ni,h,m);
277 beq(end+1,1) = 0;
278 end
279 end, end, end, end, end, end, end
280
281 % 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];
282 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
283 Aeq(end+1,:) = deltap2(i,ni,h,j,nj,k,m) - deltap2(j,nj,k,i,ni,h,m);
284 beq(end+1,1) = 0;
285 end, end, end, end, end, end, end
286
287 % 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];
288 for j = 1:M, for k =1:K(j), for nj = 1+(0:N), for i = 1:M, for m = 1:MR
289 if i~=j
290 % LHS
291 Aeq(end+1,:) = deltap2(j,nj,k,j,nj,k,m);
292 beq(end+1,1) = 0;
293 % RHS
294 for ni = 1+(0:(N-nj)) %% added +1
295 for h = 1:K(i)
296 Aeq(end,:) = Aeq(end,:) - deltap2(j,nj,k,i,ni,h,m);
297 end
298 end
299 end
300 end, end, end, end, end
301
302 %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];
303 for j = 1:M, for i = 1:M, for ki = 1:K(i)
304 % LHS
305 Aeq(end+1,:) = deltae(i,ki);
306 % RHS
307 for nj = 1+(0:N), for kj = 1:K(j), for m = 1:MR, for ni = 1+(1:N)
308 if BB(m,i)==0
309 Aeq(end,:) = Aeq(end,:) - deltap2(j,nj,kj,i,ni,ki,m);
310 end
311 end, end, end, end
312 beq(end+1,1) = 0;
313 end, end, end
314
315 %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];
316 for i = 1:M, for k =1:K(i)
317 Aeq(end+1,:) = 0;
318 beq(end+1,1) = 0;
319 % LHS
320 for j = 1:M, for h = 1:K(i)
321 Aeq(end,:) = Aeq(end,:) + q(i,j,k,h)*deltae(i,k);
322 end, end
323 % RHS
324 for j = 1:M, for h = 1:K(i)
325 Aeq(end,:) = Aeq(end,:) - q(i,j,h,k)*deltae(i,h);
326 end, end
327 end, end
328
329 %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];
330 for j = 1:M, for k =1:K(j), for nj = 1+(0:F(j)), for m = 1:MR
331 Aeq(end+1,:) = 0;
332 beq(end+1,1) = 0;
333 % LHS
334 for i = 1:M, for ni = 1+(1:F(i)), for ki = 1:K(i)
335 Aeq(end,:) = Aeq(end,:) + (ni-1)*deltap2(j,nj,k,i,ni,ki,m); % recaled back ni
336 end, end, end
337 % RHS
338 Aeq(end,:) = Aeq(end,:) - N*deltap2(j,nj,k,j,nj,k,m);
339 end, end, end, end
340
341 %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;
342 Aeq(end+1,:) = 0;
343 % LHS
344 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)
345 Aeq(end,:) = Aeq(end,:) + (ni-1)*(nj-1)*deltap2(j,nj,kj,i,ni,ki,m); % rescaled back ni and nj
346 end, end, end, end, end, end, end
347 % RHS
348 beq(end+1,1) = N^2;
349
350 %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}
351 % 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]
352 % = 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]
353 % + 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]
354 % + 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] ;
355
356 f = -1; % never equal to FC station
357 for i = 1:M, for u = 1:K(i)
358 if i~=f
359 Aeq(end+1,:)=0;
360 beq(end+1,1) = 0;
361 % LHS
362 for j = 1:M, for nj = 1+(1:F(j)), for k = 1:K(j), for h = 1:K(j), for m = 1:MR
363 if j~=i && j~=f && BB(m,j)==0
364 Aeq(end,:)= Aeq(end,:) + q(j,i,k,h)*deltap2(j,nj,k, i,1+0,u, m);
365 end
366 end, end, end, end, end
367
368 for j = 1:M, for nj = 1+(1:F(j)), for k = 1:K(j), for h = 1:K(j), for m = 1:MR
369 if j~=i && j==f && MM(m,1)~=i
370 Aeq(end,:)= Aeq(end,:) + q(j,i,k,h)*deltap2(j,nj,k, i,1+0,u, m);
371 end
372 end, end, end, end, end
373 % RHS
374 for j = 1:M, for nj = 1+(0:F(j)), for k = 1:K(i), for h = 1:K(j), for m = 1:MR
375 if j~=i && j~=f && BB(m,i)==0
376 Aeq(end,:)= Aeq(end,:) - q(i,j,k,u)*deltap2(j,nj,h, i,1+0,k, m);
377 end
378 end, end, end, end, end
379
380 for j = 1:M, for nj = 1+(0:F(j)), for k = 1:K(i), for h = 1:K(j), for m = 1:MR
381 if j~=i && j==f && BB(m,i)==0 && (nj-1)<F(j) % rescaled back nj
382 Aeq(end,:)= Aeq(end,:) - q(i,j,k,u)*deltap2(j,nj,h, i,1+1,k, m);
383 end
384 end, end, end, end, end
385
386 for j = 1:M, for nj = 1+(0:F(j)), for y = 1:K(j), for m = 1:MR
387 if j~=i && j==f && BB(m,i)==1 && (nj-1)==F(j) && MM(m,1)==i % rescaled back nj
388 for p = 1:K(f), for w = 1:M
389 if w~=f && w~=i
390 Aeq(end,:)= Aeq(end,:) - q(f,w,y,p)*deltap2(f,nj,y, i,1+1,u, m);
391 end
392 end, end
393 end
394 end, end, end, end
395 end % if
396 end, end
397
398 % subject to THM3 {i in 1..M, ni in 0..(F[i]-1): i<>f}:
399 % 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]
400 % + 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]
401 % = 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]
402 % + 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]
403 % + 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] ;
404 for i = 1:M, for ni = 1+(0:(F(i)-1))
405 if i~=f
406 Aeq(end+1,:)=0;
407 beq(end+1,1) = 0;
408 % LHS
409 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
410 if j~=i && j~=f && BB(m,j)==0
411 Aeq(end,:)= Aeq(end,:) + q(j,i,k,h)*deltap2(j,nj,k, i,ni,u, m);
412 end
413 end, end, end, end, end, end
414 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
415 if j~=i && j==f && MM(m,1)~=i
416 Aeq(end,:)= Aeq(end,:) + q(j,i,k,h)*deltap2(j,nj,k, i,ni,u, m);
417 end
418 end, end, end, end, end, end
419 % RHS
420 for j = 1:M, for nj = 1+(0:F(j)), for k = 1:K(i), for u = 1:K(j), for m = 1:MR
421 if j~=i && j~=f && BB(m,i)==0
422 for h = 1:K(i)
423 Aeq(end,:)= Aeq(end,:) - q(i,j,k,h)*deltap2(j,nj,u, i,ni+1,k, m);
424 end
425 end
426 end, end, end, end, end
427 for j = 1:M, for nj = 1+(0:F(j)), for k = 1:K(i), for u = 1:K(j), for m = 1:MR
428 if j~=i && j==f && BB(m,i)==0 && (nj-1)<F(j) % rescaled back nj
429 for h = 1:K(i)
430 Aeq(end,:)= Aeq(end,:) - q(i,j,k,h)*deltap2(j,nj,u, i,ni+1,k, m);
431 end
432 end
433 end, end, end, end, end
434 for j = 1:M, for nj = 1+(0:F(j)), for k = 1:K(i), for u = 1:K(j), for m = 1:MR
435 if j~=i && j==f && BB(m,i)==1 && (nj-1)==F(j) && MM(m,1)==i % rescaled back nj
436 for p = 1:K(f), for w = 1:M
437 if w~=f && w~=i
438 Aeq(end,:)= Aeq(end,:) - q(f,w,u,p)*deltap2(f,nj,u, i,ni+1,k, m);
439 end
440 end, end
441 end
442 end, end, end, end, end
443 end
444 end, end
445
446 % 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]
447 % + 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]
448 % = 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];
449 for i = 1:M, for ni = 1+(0:(F(i)-1))
450 if i==f
451 Aeq(end+1,:) = 0;
452 beq(end+1,1) = 0;
453 % LHS
454 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
455 if j~=i && j~=f && BB(m,j)==0 && (ni-1) < F(i) % rescaled back ni
456 Aeq(end,:) = Aeq(end,:) + q(j,i,k,h)*deltap2(j,nj,k, i,ni,u, m);
457 end
458 end, end, end, end, end, end
459 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
460 if j~=i && j==f && (ni-1)==F(i) && MM(m,1)==j % rescaled back ni
461 for w = 1:M
462 if w~=f
463 Aeq(end,:) = Aeq(end,:) + q(j,i,k,h)*deltap2(j,nj,k, i,ni,u, m);
464 end
465 end
466 end
467 end, end, end, end, end, end
468 % RHS
469 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)
470 if j~=i && (ni-1) < F(i) % rescaled back ni
471 Aeq(end,:) = Aeq(end,:) - q(i,j,k,h)*deltap2(j,nj,u, i,ni+1,k, 1);
472 end
473 end, end, end, end, end
474 end % if
475 end, end
476
477 % 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]
478 % >= 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]);
479 for j = 1:M, for k = 1:K(j), for i = 1:M, for m = 1:MR
480 A(end+1,:) = 0; % <= inequality
481 b(end+1,1) = 0;
482 % LHS with sign swapped since >= in GLPK
483 for t = 1:M, for h = 1:K(t), for nj = 1+(0:N), for nt = 1+(0:N)
484 A(end,:) = A(end,:) - (nt-1)*deltap2(j,nj,k,t,nt,h,m); % rescaled back nt
485 end, end, end, end
486 % RHS with sign swapped since >= in GLPK
487 for h = 1:K(i), for nj = 1+(0:N), for ni = 1+(1:N)
488 A(end,:) = A(end,:) + N*(deltap2(j,nj,k,i,ni,h,m));
489 end, end, end
490 end, end, end, end
491 end
492end
493
Generated by doxygen 1.9.8