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