LINE Solver
MATLAB API documentation
Loading...
Searching...
No Matches
qrf_bas.m
1function [result, x, fval, exitflag] = qrf_bas(params, objective, sense)
2% QRF_BAS - Quadratic Reduction Framework for BAS (Blocking-After-Service) networks
3%
4% MATLAB port of the AMPL model qrboundsbas_skel.mod
5%
6% Usage:
7% [result, x, fval, exitflag] = qrf_bas(params)
8% [result, x, fval, exitflag] = qrf_bas(params, objective)
9% [result, x, fval, exitflag] = qrf_bas(params, objective, sense)
10%
11% Inputs:
12% params - Structure with model parameters:
13% .M - Number of queues
14% .N - Total population
15% .f - Index of finite capacity queue (1-based)
16% .F - [M x 1] Capacity of each queue
17% .K - [M x 1] Number of phases for each queue
18% .mu - {M x 1} cell, each mu{i} is K(i) x K(i) completion rates
19% .v - {M x 1} cell, each v{i} is K(i) x K(i) background rates
20% .r - [M x M] Routing probabilities
21% .MR - Number of blocking configurations
22% .BB - [MR x M] Blocking state (0/1)
23% .MM - [MR x 2] Blocking order (queue indices)
24% .ZZ - [MR x 1] Number of blocked queues in each config
25% .ZM - Maximum blocking depth
26% .MM1 - [MR x M] Extended blocking order info
27%
28% .lpAlgorithm - (optional) linprog Algorithm string;
29% defaults to 'interior-point-legacy' (see note at the
30% LP solve for why, and for the R2025a HiGHS caveat)
31%
32% objective - (optional) 'U1min' (default), 'U1max', or queue index 1..M
33% sense - (optional) 'min' (default) or 'max'
34%
35% Outputs:
36% result - Structure with results:
37% .U - [M x 1] Utilization of each queue
38% .e - [M x max(K)] Effective utilization by phase
39% x - Raw solution vector
40% fval - Objective function value
41% exitflag - Solver exit flag
42
43 if nargin < 2 || isempty(objective)
44 objective = 'U1min';
45 end
46 if nargin < 3 || isempty(sense)
47 sense = 'min';
48 end
49
50 % Extract parameters
51 M = params.M;
52 N = params.N;
53 f = params.f; % finite capacity queue index
54 F = params.F(:);
55 K = params.K(:);
56 mu = params.mu;
57 v = params.v;
58 r = params.r;
59 MR = params.MR;
60 BB = params.BB;
61 MM = params.MM;
62 ZZ = params.ZZ(:);
63 ZM = params.ZM;
64 MM1 = params.MM1;
65 if isfield(params, 'verbose')
66 verbose = params.verbose;
67 else
68 verbose = true;
69 end
70
71 % Compute transition rates q(i,j,k,h)
72 % q{i,j} is a K(i) x K(i) array (not load-dependent for BAS)
73 q = cell(M, M);
74 for i = 1:M
75 for j = 1:M
76 q{i,j} = zeros(K(i), K(i));
77 for ki = 1:K(i)
78 for hi = 1:K(i)
79 if j ~= i
80 q{i,j}(ki, hi) = r(i,j) * mu{i}(ki, hi);
81 else
82 q{i,j}(ki, hi) = v{i}(ki, hi) + r(i,i) * mu{i}(ki, hi);
83 end
84 end
85 end
86 end
87 end
88
89 %% Build variable indexing
90 % p2(j, nj, kj, i, ni, hi, m) for j in 1:M, nj in 0:N, kj in 1:K(j),
91 % i in 1:M, ni in 0:N, hi in 1:K(i), m in 1:MR
92 % e(i, ki) for i in 1:M, ki in 1:K(i)
93
94 if verbose; fprintf('Building variable index map...\n'); end
95 varCount = 0;
96
97 % p2 variables
98 p2idx = cell(M, 1);
99 for j = 1:M
100 p2idx{j} = cell(N+1, K(j), M, MR);
101 for nj = 0:N
102 for kj = 1:K(j)
103 for i = 1:M
104 for m = 1:MR
105 p2idx{j}{nj+1, kj, i, m} = zeros(N+1, K(i));
106 for ni = 0:N
107 for hi = 1:K(i)
108 varCount = varCount + 1;
109 p2idx{j}{nj+1, kj, i, m}(ni+1, hi) = varCount;
110 end
111 end
112 end
113 end
114 end
115 end
116 end
117
118 % e variables
119 eidx = cell(M, 1);
120 for i = 1:M
121 eidx{i} = zeros(K(i), 1);
122 for ki = 1:K(i)
123 varCount = varCount + 1;
124 eidx{i}(ki) = varCount;
125 end
126 end
127
128 nVars = varCount;
129 if verbose; fprintf('Total variables: %d\n', nVars); end
130
131 %% Build constraints
132 % Sparse triplet accumulation. The previous builder used
133 % row = zeros(1,nVars); row(idx) = row(idx) + v; Aeq = [Aeq; row];
134 % which is O(nrows^2) in memory traffic: every append reallocated and copied
135 % the whole matrix, and every row allocated nVars doubles. For the BAS
136 % example (nVars ~6e4, ~1.1e5 rows) that does not terminate in practice.
137 % We emit (row,col,val) triplets and build one sparse matrix at the end.
138 % sparse() SUMS duplicate (i,j) entries, which reproduces the
139 % row(idx)=row(idx)+v accumulation exactly. The plain assignments below
140 % (SYMMETRY +-1, UEFF row(idx_e)=-1, THM2 row(idx_diag)=-N, ONE/MARGINALS
141 % =1) are each the FIRST write to their index within their row, so emitting
142 % them as triplets is equivalent; SYMMETRY additionally guards idx1~=idx2.
143 eqNnzCap = 1048576;
144 eqI = zeros(eqNnzCap, 1); eqJ = zeros(eqNnzCap, 1); eqV = zeros(eqNnzCap, 1);
145 eqNnz = 0; nEq = 0;
146 inNnzCap = 65536;
147 inI = zeros(inNnzCap, 1); inJ = zeros(inNnzCap, 1); inV = zeros(inNnzCap, 1);
148 inNnz = 0; nIn = 0;
149 beq = zeros(65536, 1);
150 bineq = zeros(4096, 1);
151 rCap = 4096; rI = zeros(rCap, 1); rV = zeros(rCap, 1); rN = 0;
152
153 % Helper functions
154 getP2Idx = @(j, nj, kj, i, ni, hi, m) p2idx{j}{nj+1, kj, i, m}(ni+1, hi);
155 getEIdx = @(i, ki) eidx{i}(ki);
156
157 if verbose; fprintf('Building constraints...\n'); end
158
159 %% Initialize bounds
160 lb = zeros(nVars, 1);
161 ub = inf(nVars, 1);
162
163 %% ZERO constraints - fix infeasible states
164 if verbose; fprintf(' ZERO constraints...\n'); end
165 for j = 1:M
166 for nj = 0:N
167 for kj = 1:K(j)
168 for i = 1:M
169 for ni = 0:N
170 for hi = 1:K(i)
171 for m = 1:MR
172 idx = getP2Idx(j, nj, kj, i, ni, hi, m);
173
174 % ZERO1: i==j, nj==ni, h<>k
175 if i == j && nj == ni && hi ~= kj
176 ub(idx) = 0;
177 end
178
179 % ZERO2: i==j, nj<>ni
180 if i == j && nj ~= ni
181 ub(idx) = 0;
182 end
183
184 % ZERO3: i<>j, nj+ni > N
185 if i ~= j && nj + ni > N
186 ub(idx) = 0;
187 end
188
189 % ZERO6: nj > F(j)
190 if nj > F(j)
191 ub(idx) = 0;
192 end
193
194 % ZERO5: BB(m,j)==1 and nj==0
195 if m >= 2 && BB(m, j) == 1 && nj == 0
196 ub(idx) = 0;
197 end
198
199 % ZERO7: BB(m,j)==1 and i<>j and i<>f and ni+nj+F(f)>N
200 if m >= 2 && BB(m, j) == 1 && i ~= j && i ~= f && ni + nj + F(f) > N
201 ub(idx) = 0;
202 end
203
204 % ZERO8: finite queue not at capacity in blocking config
205 if j == f && nj >= 1 && nj <= F(f)-1 && m >= 2
206 ub(idx) = 0;
207 end
208 end
209 end
210 end
211 end
212 end
213 end
214 end
215
216 % ZERO4: For m>=2 and j<>f, p2(j,nj,k,f,nf,h,m)=0 when nf < F(f)
217 for j = 1:M
218 if j == f
219 continue;
220 end
221 for nj = 0:N
222 for kj = 1:K(j)
223 for m = 2:MR
224 for nf = 0:(F(f)-1)
225 for hf = 1:K(f)
226 idx = getP2Idx(j, nj, kj, f, nf, hf, m);
227 ub(idx) = 0;
228 end
229 end
230 end
231 end
232 end
233 end
234
235 %% ONE: Normalization
236 if verbose; fprintf(' ONE constraints...\n'); end
237 for j = 1:M
238 rN = 0;
239 for nj = 0:N
240 for kj = 1:K(j)
241 for m = 1:MR
242 idx = getP2Idx(j, nj, kj, j, nj, kj, m);
243 rN=rN+1; if rN>rCap, rCap=2*rCap; rI(rCap)=0; rV(rCap)=0; end; rI(rN)=idx; rV(rN)=1;
244 end
245 end
246 end
247 nEq=nEq+1; if eqNnz+rN>numel(eqI), eqI(2*(eqNnz+rN))=0; eqJ(2*(eqNnz+rN))=0; eqV(2*(eqNnz+rN))=0; end; eqI(eqNnz+(1:rN))=nEq; eqJ(eqNnz+(1:rN))=rI(1:rN); eqV(eqNnz+(1:rN))=rV(1:rN); eqNnz=eqNnz+rN; if nEq>numel(beq), beq(2*nEq)=0; end
248 beq(nEq)=1;
249 end
250
251 %% SYMMETRY
252 if verbose; fprintf(' SYMMETRY constraints...\n'); end
253 for j = 1:M
254 for nj = 0:min(N, F(j))
255 for kj = 1:K(j)
256 for i = 1:M
257 if i <= j
258 continue;
259 end
260 for ni = 0:min(N, F(i))
261 if i ~= j && nj + ni > N
262 continue;
263 end
264 for hi = 1:K(i)
265 for m = 1:MR
266 idx1 = getP2Idx(j, nj, kj, i, ni, hi, m);
267 idx2 = getP2Idx(i, ni, hi, j, nj, kj, m);
268 if ub(idx1) == 0 && ub(idx2) == 0
269 continue;
270 end
271 if idx1 ~= idx2
272 rN = 0;
273 rN=rN+1; if rN>rCap, rCap=2*rCap; rI(rCap)=0; rV(rCap)=0; end; rI(rN)=idx1; rV(rN)=1;
274 rN=rN+1; if rN>rCap, rCap=2*rCap; rI(rCap)=0; rV(rCap)=0; end; rI(rN)=idx2; rV(rN)=-1;
275 nEq=nEq+1; if eqNnz+rN>numel(eqI), eqI(2*(eqNnz+rN))=0; eqJ(2*(eqNnz+rN))=0; eqV(2*(eqNnz+rN))=0; end; eqI(eqNnz+(1:rN))=nEq; eqJ(eqNnz+(1:rN))=rI(1:rN); eqV(eqNnz+(1:rN))=rV(1:rN); eqNnz=eqNnz+rN; if nEq>numel(beq), beq(2*nEq)=0; end
276 beq(nEq)=0;
277 end
278 end
279 end
280 end
281 end
282 end
283 end
284 end
285
286 %% MARGINALS
287 if verbose; fprintf(' MARGINALS constraints...\n'); end
288 for j = 1:M
289 for kj = 1:K(j)
290 for nj = 0:min(N, F(j))
291 for i = 1:M
292 if i == j
293 continue;
294 end
295 for m = 1:MR
296 rN = 0;
297 idx_diag = getP2Idx(j, nj, kj, j, nj, kj, m);
298 rN=rN+1; if rN>rCap, rCap=2*rCap; rI(rCap)=0; rV(rCap)=0; end; rI(rN)=idx_diag; rV(rN)=1;
299 for ni = 0:min(N-nj, F(i))
300 for hi = 1:K(i)
301 idx = getP2Idx(j, nj, kj, i, ni, hi, m);
302 rN=rN+1; if rN>rCap, rCap=2*rCap; rI(rCap)=0; rV(rCap)=0; end; rI(rN)=idx; rV(rN)=-(1);
303 end
304 end
305 nEq=nEq+1; if eqNnz+rN>numel(eqI), eqI(2*(eqNnz+rN))=0; eqJ(2*(eqNnz+rN))=0; eqV(2*(eqNnz+rN))=0; end; eqI(eqNnz+(1:rN))=nEq; eqJ(eqNnz+(1:rN))=rI(1:rN); eqV(eqNnz+(1:rN))=rV(1:rN); eqNnz=eqNnz+rN; if nEq>numel(beq), beq(2*nEq)=0; end
306 beq(nEq)=0;
307 end
308 end
309 end
310 end
311 end
312
313 %% UEFF: e(i,ki) = sum of p2 where queue i is not blocked
314 if verbose; fprintf(' UEFF constraints...\n'); end
315 for i = 1:M
316 for ki = 1:K(i)
317 rN = 0;
318 idx_e = getEIdx(i, ki);
319 rN=rN+1; if rN>rCap, rCap=2*rCap; rI(rCap)=0; rV(rCap)=0; end; rI(rN)=idx_e; rV(rN)=-1;
320 for j = 1:M
321 for nj = 0:min(N, F(j))
322 for kj = 1:K(j)
323 for m = 1:MR
324 if BB(m, i) == 0
325 for ni = 1:min(N, F(i))
326 idx = getP2Idx(j, nj, kj, i, ni, ki, m);
327 rN=rN+1; if rN>rCap, rCap=2*rCap; rI(rCap)=0; rV(rCap)=0; end; rI(rN)=idx; rV(rN)=1;
328 end
329 end
330 end
331 end
332 end
333 end
334 nEq=nEq+1; if eqNnz+rN>numel(eqI), eqI(2*(eqNnz+rN))=0; eqJ(2*(eqNnz+rN))=0; eqV(2*(eqNnz+rN))=0; end; eqI(eqNnz+(1:rN))=nEq; eqJ(eqNnz+(1:rN))=rI(1:rN); eqV(eqNnz+(1:rN))=rV(1:rN); eqNnz=eqNnz+rN; if nEq>numel(beq), beq(2*nEq)=0; end
335 beq(nEq)=0;
336 end
337 end
338
339 %% THM1: Phase balance (Theorem 1)
340 % sum {j, h: j<>i or h<>k} q(i,j,k,h)*e(i,k) = sum {j, h: j<>i or h<>k} q(i,j,h,k)*e(i,h)
341 if verbose; fprintf(' THM1 (Phase balance) constraints...\n'); end
342 for i = 1:M
343 for ki = 1:K(i)
344 rN = 0;
345 % LHS
346 for j = 1:M
347 for hi = 1:K(i)
348 if j ~= i || hi ~= ki
349 idx_e = getEIdx(i, ki);
350 rN=rN+1; if rN>rCap, rCap=2*rCap; rI(rCap)=0; rV(rCap)=0; end; rI(rN)=idx_e; rV(rN)=q{i,j}(ki, hi);
351 end
352 end
353 end
354 % RHS (subtract)
355 for j = 1:M
356 for hi = 1:K(i)
357 if j ~= i || hi ~= ki
358 idx_e = getEIdx(i, hi);
359 rN=rN+1; if rN>rCap, rCap=2*rCap; rI(rCap)=0; rV(rCap)=0; end; rI(rN)=idx_e; rV(rN)=-(q{i,j}(hi, ki));
360 end
361 end
362 end
363 nEq=nEq+1; if eqNnz+rN>numel(eqI), eqI(2*(eqNnz+rN))=0; eqJ(2*(eqNnz+rN))=0; eqV(2*(eqNnz+rN))=0; end; eqI(eqNnz+(1:rN))=nEq; eqJ(eqNnz+(1:rN))=rI(1:rN); eqV(eqNnz+(1:rN))=rV(1:rN); eqNnz=eqNnz+rN; if nEq>numel(beq), beq(2*nEq)=0; end
364 beq(nEq)=0;
365 end
366 end
367
368 %% THM2: Population constraint (Theorem 2)
369 if verbose; fprintf(' THM2 (Population) constraints...\n'); end
370 for j = 1:M
371 for kj = 1:K(j)
372 for nj = 0:F(j)
373 for m = 1:MR
374 rN = 0;
375 % RHS: -N * p2(j,nj,kj,j,nj,kj,m)
376 idx_diag = getP2Idx(j, nj, kj, j, nj, kj, m);
377 rN=rN+1; if rN>rCap, rCap=2*rCap; rI(rCap)=0; rV(rCap)=0; end; rI(rN)=idx_diag; rV(rN)=-N;
378 % LHS: sum
379 for i = 1:M
380 for ni = 1:F(i)
381 for ki = 1:K(i)
382 idx = getP2Idx(j, nj, kj, i, ni, ki, m);
383 rN=rN+1; if rN>rCap, rCap=2*rCap; rI(rCap)=0; rV(rCap)=0; end; rI(rN)=idx; rV(rN)=ni;
384 end
385 end
386 end
387 nEq=nEq+1; if eqNnz+rN>numel(eqI), eqI(2*(eqNnz+rN))=0; eqJ(2*(eqNnz+rN))=0; eqV(2*(eqNnz+rN))=0; end; eqI(eqNnz+(1:rN))=nEq; eqJ(eqNnz+(1:rN))=rI(1:rN); eqV(eqNnz+(1:rN))=rV(1:rN); eqNnz=eqNnz+rN; if nEq>numel(beq), beq(2*nEq)=0; end
388 beq(nEq)=0;
389 end
390 end
391 end
392 end
393
394 %% COR1: Second moment constraint (Corollary to Theorem 2)
395 % sum_{m,i,j,nj,ni,ki,kj} ni*nj*p2(j,nj,kj,i,ni,ki,m) = N^2
396 if verbose; fprintf(' COR1 (Second moment) constraint...\n'); end
397 rN = 0;
398 for m = 1:MR
399 for i = 1:M
400 for j = 1:M
401 for nj = 1:F(j)
402 for ni = 1:F(i)
403 for ki = 1:K(i)
404 for kj = 1:K(j)
405 idx = getP2Idx(j, nj, kj, i, ni, ki, m);
406 rN=rN+1; if rN>rCap, rCap=2*rCap; rI(rCap)=0; rV(rCap)=0; end; rI(rN)=idx; rV(rN)=ni * nj;
407 end
408 end
409 end
410 end
411 end
412 end
413 end
414 nEq=nEq+1; if eqNnz+rN>numel(eqI), eqI(2*(eqNnz+rN))=0; eqJ(2*(eqNnz+rN))=0; eqV(2*(eqNnz+rN))=0; end; eqI(eqNnz+(1:rN))=nEq; eqJ(eqNnz+(1:rN))=rI(1:rN); eqV(eqNnz+(1:rN))=rV(1:rN); eqNnz=eqNnz+rN; if nEq>numel(beq), beq(2*nEq)=0; end
415 beq(nEq)=N^2;
416
417 %% THM30: Marginal balance for ni=0 (per phase), i<>f
418 if verbose; fprintf(' THM30 (Marginal balance ni=0) constraints...\n'); end
419 for i = 1:M
420 if i == f
421 continue;
422 end
423 for ui = 1:K(i)
424 rN = 0;
425 % LHS: arrivals from j<>i,j<>f with BB(m,j)==0
426 for j = 1:M
427 if j == i || j == f
428 continue;
429 end
430 for nj = 1:F(j)
431 for kj = 1:K(j)
432 for hj = 1:K(j)
433 for m = 1:MR
434 if BB(m, j) == 0
435 idx = getP2Idx(j, nj, kj, i, 0, ui, m);
436 rN=rN+1; if rN>rCap, rCap=2*rCap; rI(rCap)=0; rV(rCap)=0; end; rI(rN)=idx; rV(rN)=q{j,i}(kj, hj);
437 end
438 end
439 end
440 end
441 end
442 end
443 % LHS: arrivals from j==f with MM(m,1)<>i
444 for nj = 1:F(f)
445 for kj = 1:K(f)
446 for hj = 1:K(f)
447 for m = 1:MR
448 if MM(m, 1) ~= i
449 idx = getP2Idx(f, nj, kj, i, 0, ui, m);
450 rN=rN+1; if rN>rCap, rCap=2*rCap; rI(rCap)=0; rV(rCap)=0; end; rI(rN)=idx; rV(rN)=q{f,i}(kj, hj);
451 end
452 end
453 end
454 end
455 end
456
457 % RHS: departures from i at ni=1 to j<>i,j<>f with BB(m,i)==0
458 for j = 1:M
459 if j == i || j == f
460 continue;
461 end
462 for nj = 0:F(j)
463 for ki = 1:K(i)
464 for hj = 1:K(j)
465 for m = 1:MR
466 if BB(m, i) == 0
467 idx = getP2Idx(j, nj, hj, i, 1, ki, m);
468 rN=rN+1; if rN>rCap, rCap=2*rCap; rI(rCap)=0; rV(rCap)=0; end; rI(rN)=idx; rV(rN)=-(q{i,j}(ki, ui));
469 end
470 end
471 end
472 end
473 end
474 end
475 % RHS: departures to j==f with BB(m,i)==0 and nj<F(f)
476 for nj = 0:(F(f)-1)
477 for ki = 1:K(i)
478 for hj = 1:K(f)
479 for m = 1:MR
480 if BB(m, i) == 0
481 idx = getP2Idx(f, nj, hj, i, 1, ki, m);
482 rN=rN+1; if rN>rCap, rCap=2*rCap; rI(rCap)=0; rV(rCap)=0; end; rI(rN)=idx; rV(rN)=-(q{i,f}(ki, ui));
483 end
484 end
485 end
486 end
487 end
488 % RHS: unblocking when BB(m,i)==1 and MM(m,1)==i
489 for m = 1:MR
490 if BB(m, i) == 1 && MM(m, 1) == i
491 for kf = 1:K(f)
492 for pf = 1:K(f)
493 for w = 1:M
494 if w ~= f && w ~= i
495 idx = getP2Idx(f, F(f), kf, i, 1, ui, m);
496 rN=rN+1; if rN>rCap, rCap=2*rCap; rI(rCap)=0; rV(rCap)=0; end; rI(rN)=idx; rV(rN)=-(q{f,w}(kf, pf));
497 end
498 end
499 end
500 end
501 end
502 end
503
504 nEq=nEq+1; if eqNnz+rN>numel(eqI), eqI(2*(eqNnz+rN))=0; eqJ(2*(eqNnz+rN))=0; eqV(2*(eqNnz+rN))=0; end; eqI(eqNnz+(1:rN))=nEq; eqJ(eqNnz+(1:rN))=rI(1:rN); eqV(eqNnz+(1:rN))=rV(1:rN); eqNnz=eqNnz+rN; if nEq>numel(beq), beq(2*nEq)=0; end
505 beq(nEq)=0;
506 end
507 end
508
509 %% THM3: Marginal balance for ni in 1:F(i)-1, i<>f
510 if verbose; fprintf(' THM3 (Marginal balance) constraints...\n'); end
511 for i = 1:M
512 if i == f
513 continue;
514 end
515 for ni = 1:(F(i)-1)
516 rN = 0;
517 % LHS: arrivals from j<>i,j<>f with BB(m,j)==0
518 for j = 1:M
519 if j == i || j == f
520 continue;
521 end
522 for nj = 1:F(j)
523 for kj = 1:K(j)
524 for hj = 1:K(j)
525 for ui = 1:K(i)
526 for m = 1:MR
527 if BB(m, j) == 0
528 idx = getP2Idx(j, nj, kj, i, ni, ui, m);
529 rN=rN+1; if rN>rCap, rCap=2*rCap; rI(rCap)=0; rV(rCap)=0; end; rI(rN)=idx; rV(rN)=q{j,i}(kj, hj);
530 end
531 end
532 end
533 end
534 end
535 end
536 end
537 % LHS: arrivals from j==f with MM(m,1)<>i
538 for nj = 1:F(f)
539 for kj = 1:K(f)
540 for hj = 1:K(f)
541 for ui = 1:K(i)
542 for m = 1:MR
543 if MM(m, 1) ~= i
544 idx = getP2Idx(f, nj, kj, i, ni, ui, m);
545 rN=rN+1; if rN>rCap, rCap=2*rCap; rI(rCap)=0; rV(rCap)=0; end; rI(rN)=idx; rV(rN)=q{f,i}(kj, hj);
546 end
547 end
548 end
549 end
550 end
551 end
552
553 % RHS: departures from i at ni+1 to j<>i,j<>f with BB(m,i)==0
554 for j = 1:M
555 if j == i || j == f
556 continue;
557 end
558 for nj = 0:F(j)
559 for ki = 1:K(i)
560 for hi = 1:K(i)
561 for uj = 1:K(j)
562 for m = 1:MR
563 if BB(m, i) == 0
564 idx = getP2Idx(j, nj, uj, i, ni+1, ki, m);
565 rN=rN+1; if rN>rCap, rCap=2*rCap; rI(rCap)=0; rV(rCap)=0; end; rI(rN)=idx; rV(rN)=-(q{i,j}(ki, hi));
566 end
567 end
568 end
569 end
570 end
571 end
572 end
573 % RHS: departures to j==f with BB(m,i)==0 and nj<F(f)
574 for nj = 0:(F(f)-1)
575 for ki = 1:K(i)
576 for hi = 1:K(i)
577 for uj = 1:K(f)
578 for m = 1:MR
579 if BB(m, i) == 0
580 idx = getP2Idx(f, nj, uj, i, ni+1, ki, m);
581 rN=rN+1; if rN>rCap, rCap=2*rCap; rI(rCap)=0; rV(rCap)=0; end; rI(rN)=idx; rV(rN)=-(q{i,f}(ki, hi));
582 end
583 end
584 end
585 end
586 end
587 end
588 % RHS: unblocking when BB(m,i)==1 and MM(m,1)==i
589 for m = 1:MR
590 if BB(m, i) == 1 && MM(m, 1) == i
591 for ki = 1:K(i)
592 for kf = 1:K(f)
593 for pf = 1:K(f)
594 for w = 1:M
595 if w ~= f && w ~= i
596 idx = getP2Idx(f, F(f), kf, i, ni+1, ki, m);
597 rN=rN+1; if rN>rCap, rCap=2*rCap; rI(rCap)=0; rV(rCap)=0; end; rI(rN)=idx; rV(rN)=-(q{f,w}(kf, pf));
598 end
599 end
600 end
601 end
602 end
603 end
604 end
605
606 nEq=nEq+1; if eqNnz+rN>numel(eqI), eqI(2*(eqNnz+rN))=0; eqJ(2*(eqNnz+rN))=0; eqV(2*(eqNnz+rN))=0; end; eqI(eqNnz+(1:rN))=nEq; eqJ(eqNnz+(1:rN))=rI(1:rN); eqV(eqNnz+(1:rN))=rV(1:rN); eqNnz=eqNnz+rN; if nEq>numel(beq), beq(2*nEq)=0; end
607 beq(nEq)=0;
608 end
609 end
610
611 %% THM3f: Marginal balance for i==f
612 if verbose; fprintf(' THM3f (Marginal balance for finite queue) constraints...\n'); end
613 for ni = 0:(F(f)-1)
614 rN = 0;
615 % LHS: arrivals from j<>f with BB(m,j)==0 and ni<F(f)
616 if ni < F(f)
617 for j = 1:M
618 if j == f
619 continue;
620 end
621 for nj = 1:F(j)
622 for kj = 1:K(j)
623 for hj = 1:K(j)
624 for uf = 1:K(f)
625 for m = 1:MR
626 if BB(m, j) == 0
627 idx = getP2Idx(j, nj, kj, f, ni, uf, m);
628 rN=rN+1; if rN>rCap, rCap=2*rCap; rI(rCap)=0; rV(rCap)=0; end; rI(rN)=idx; rV(rN)=q{j,f}(kj, hj);
629 end
630 end
631 end
632 end
633 end
634 end
635 end
636 end
637
638 % RHS: departures from f at ni+1 to j<>f (only m=1, no blocking)
639 for j = 1:M
640 if j == f
641 continue;
642 end
643 for nj = 0:F(j)
644 for kf = 1:K(f)
645 for hf = 1:K(f)
646 for uj = 1:K(j)
647 if ni < F(f)
648 idx = getP2Idx(j, nj, uj, f, ni+1, kf, 1);
649 rN=rN+1; if rN>rCap, rCap=2*rCap; rI(rCap)=0; rV(rCap)=0; end; rI(rN)=idx; rV(rN)=-(q{f,j}(kf, hf));
650 end
651 end
652 end
653 end
654 end
655 end
656
657 nEq=nEq+1; if eqNnz+rN>numel(eqI), eqI(2*(eqNnz+rN))=0; eqJ(2*(eqNnz+rN))=0; eqV(2*(eqNnz+rN))=0; end; eqI(eqNnz+(1:rN))=nEq; eqJ(eqNnz+(1:rN))=rI(1:rN); eqV(eqNnz+(1:rN))=rV(1:rN); eqNnz=eqNnz+rN; if nEq>numel(beq), beq(2*nEq)=0; end
658 beq(nEq)=0;
659 end
660
661 %% THM3I: Blocking depth balance (Theorem 4)
662 if verbose; fprintf(' THM3I (Blocking depth balance) constraints...\n'); end
663 for z = 0:(ZM-1)
664 rN = 0;
665 % LHS: arrivals to f at F(f) from j<>f with BB(m,j)==0 and ZZ(m)==z
666 for j = 1:M
667 if j == f
668 continue;
669 end
670 for nj = 1:F(j)
671 for kj = 1:K(j)
672 for hj = 1:K(j)
673 for uf = 1:K(f)
674 for m = 1:MR
675 if BB(m, j) == 0 && ZZ(m) == z
676 idx = getP2Idx(j, nj, kj, f, F(f), uf, m);
677 rN=rN+1; if rN>rCap, rCap=2*rCap; rI(rCap)=0; rV(rCap)=0; end; rI(rN)=idx; rV(rN)=q{j,f}(kj, hj);
678 end
679 end
680 end
681 end
682 end
683 end
684 end
685 % RHS: departures from f to j<>f with ZZ(m)==z+1
686 for j = 1:M
687 if j == f
688 continue;
689 end
690 for nj = 0:F(j)
691 for kf = 1:K(f)
692 for hf = 1:K(f)
693 for uj = 1:K(j)
694 for m = 1:MR
695 if ZZ(m) == z+1
696 idx = getP2Idx(j, nj, uj, f, F(f), kf, m);
697 rN=rN+1; if rN>rCap, rCap=2*rCap; rI(rCap)=0; rV(rCap)=0; end; rI(rN)=idx; rV(rN)=-(q{f,j}(kf, hf));
698 end
699 end
700 end
701 end
702 end
703 end
704 end
705
706 nEq=nEq+1; if eqNnz+rN>numel(eqI), eqI(2*(eqNnz+rN))=0; eqJ(2*(eqNnz+rN))=0; eqV(2*(eqNnz+rN))=0; end; eqI(eqNnz+(1:rN))=nEq; eqJ(eqNnz+(1:rN))=rI(1:rN); eqV(eqNnz+(1:rN))=rV(1:rN); eqNnz=eqNnz+rN; if nEq>numel(beq), beq(2*nEq)=0; end
707 beq(nEq)=0;
708 end
709
710 %% THM3L: Maximum blocking depth constraint
711 if verbose; fprintf(' THM3L (Max blocking depth) constraints...\n'); end
712 for m = 1:MR
713 if ZZ(m) ~= ZM - 1
714 continue;
715 end
716 rN = 0;
717 % LHS: arrivals from j<>f with BB(m,j)==0 and MM1(m,j)>0
718 for j = 1:M
719 if j == f || BB(m, j) ~= 0 || MM1(m, j) <= 0
720 continue;
721 end
722 for nj = 1:F(j)
723 for kj = 1:K(j)
724 for hj = 1:K(j)
725 for uf = 1:K(f)
726 idx = getP2Idx(j, nj, kj, f, F(f), uf, m);
727 rN=rN+1; if rN>rCap, rCap=2*rCap; rI(rCap)=0; rV(rCap)=0; end; rI(rN)=idx; rV(rN)=q{j,f}(kj, hj);
728 end
729 end
730 end
731 end
732 end
733 % RHS: uses MM1(m,j) to index into blocking configuration
734 for j = 1:M
735 if j == f || BB(m, j) ~= 0 || MM1(m, j) <= 0
736 continue;
737 end
738 mp = MM1(m, j); % blocking configuration index
739 for kf = 1:K(f)
740 for uf = 1:K(f)
741 for w = 1:M
742 if w ~= f
743 idx = getP2Idx(f, F(f), kf, f, F(f), kf, mp);
744 rN=rN+1; if rN>rCap, rCap=2*rCap; rI(rCap)=0; rV(rCap)=0; end; rI(rN)=idx; rV(rN)=-(q{f,w}(kf, uf));
745 end
746 end
747 end
748 end
749 end
750
751 nEq=nEq+1; if eqNnz+rN>numel(eqI), eqI(2*(eqNnz+rN))=0; eqJ(2*(eqNnz+rN))=0; eqV(2*(eqNnz+rN))=0; end; eqI(eqNnz+(1:rN))=nEq; eqJ(eqNnz+(1:rN))=rI(1:rN); eqV(eqNnz+(1:rN))=rV(1:rN); eqNnz=eqNnz+rN; if nEq>numel(beq), beq(2*nEq)=0; end
752 beq(nEq)=0;
753 end
754
755 %% THM4: Queue-length bound inequality (Theorem 5)
756 if verbose; fprintf(' THM4 (Queue-length bound) constraints...\n'); end
757 for j = 1:M
758 for kj = 1:K(j)
759 for i = 1:M
760 for m = 1:MR
761 rN = 0;
762 % LHS: sum_t sum_ht sum_nj sum_nt nt * p2
763 for t = 1:M
764 for ht = 1:K(t)
765 for nj = 0:F(j)
766 for nt = 1:F(t)
767 idx = getP2Idx(j, nj, kj, t, nt, ht, m);
768 rN=rN+1; if rN>rCap, rCap=2*rCap; rI(rCap)=0; rV(rCap)=0; end; rI(rN)=idx; rV(rN)=nt;
769 end
770 end
771 end
772 end
773 % RHS: -N * sum
774 for hi = 1:K(i)
775 for nj = 0:F(j)
776 for ni = 1:F(i)
777 idx = getP2Idx(j, nj, kj, i, ni, hi, m);
778 rN=rN+1; if rN>rCap, rCap=2*rCap; rI(rCap)=0; rV(rCap)=0; end; rI(rN)=idx; rV(rN)=-(N);
779 end
780 end
781 end
782 nIn=nIn+1; if inNnz+rN>numel(inI), inI(2*(inNnz+rN))=0; inJ(2*(inNnz+rN))=0; inV(2*(inNnz+rN))=0; end; inI(inNnz+(1:rN))=nIn; inJ(inNnz+(1:rN))=rI(1:rN); inV(inNnz+(1:rN))=-rV(1:rN); inNnz=inNnz+rN; if nIn>numel(bineq), bineq(2*nIn)=0; end
783 bineq(nIn)=0;
784 end
785 end
786 end
787 end
788
789 % Materialize the sparse constraint matrices from the triplets.
790 Aeq = sparse(eqI(1:eqNnz), eqJ(1:eqNnz), eqV(1:eqNnz), nEq, nVars);
791 beq = beq(1:nEq);
792 Aineq = sparse(inI(1:inNnz), inJ(1:inNnz), inV(1:inNnz), nIn, nVars);
793 bineq = bineq(1:nIn);
794
795 %% Build objective function
796 if verbose; fprintf('Building objective function...\n'); end
797 c = zeros(nVars, 1);
798
799 if ischar(objective)
800 if strcmp(objective, 'U1min') || strcmp(objective, 'U1max')
801 targetQueue = 1;
802 else
803 error('Unknown objective: %s', objective);
804 end
805 else
806 targetQueue = objective;
807 end
808
809 % Utilization = sum over m, k, n of p2(i,n,k,i,n,k,m)
810 for m = 1:MR
811 for ki = 1:K(targetQueue)
812 for ni = 1:F(targetQueue)
813 idx = getP2Idx(targetQueue, ni, ki, targetQueue, ni, ki, m);
814 c(idx) = 1;
815 end
816 end
817 end
818
819 if strcmp(sense, 'max')
820 c = -c;
821 end
822
823 %% Solve LP
824 if verbose; fprintf('Solving LP with %d variables and %d equality + %d inequality constraints...\n', ...
825 nVars, size(Aeq, 1), size(Aineq, 1)); end
826
827 % LP algorithm. R2025a's default 'dual-simplex-highs' is broken in some
828 % installs (errors "Unrecognized field name optimstatus"), so an
829 % interior-point variant is required here. On the paper's BAS instance,
830 % which is badly scaled (mu spans 1.016186 down to 2.585708e-05),
831 % 'interior-point-legacy' is markedly more accurate than 'interior-point':
832 % against the published GLPK optimum it gives |dU1min|=8.5e-07 and
833 % |dU1max|=9.2e-10, versus 4.3e-05 and 8.0e-06 for 'interior-point'. The
834 % residual is pure solver tolerance, not formulation: the minimum lands
835 % above and the maximum below the GLPK vertex, and both gaps shrink
836 % together as the solver becomes more accurate. Override via
837 % params.lpAlgorithm if a particular model needs a different method.
838 if isfield(params, 'lpAlgorithm') && ~isempty(params.lpAlgorithm)
839 lpAlgorithm = params.lpAlgorithm;
840 else
841 lpAlgorithm = 'interior-point-legacy';
842 end
843 if verbose
844 options = optimoptions('linprog', 'Display', 'final', 'Algorithm', lpAlgorithm);
845 else
846 options = optimoptions('linprog', 'Display', 'off', 'Algorithm', lpAlgorithm);
847 end
848
849 [x, fval, exitflag] = linprog(c, Aineq, bineq, Aeq, beq, lb, ub, options);
850
851 if strcmp(sense, 'max')
852 fval = -fval;
853 end
854
855 %% Extract results
856 result = struct();
857 result.objective = fval;
858 result.exitflag = exitflag;
859
860 % Whether the metric fields below can be filled is a property of the
861 % solution vector, not of exitflag. exitflag 0 (iteration limit) still
862 % returns a usable interior point -- that is the normal outcome on badly
863 % scaled instances such as the reference BAS network -- whereas a solve that
864 % breaks down returns an empty or non-finite x, which must never be
865 % consumed. Predicate on x accordingly, and say so rather than skipping
866 % silently. (Observed here: exitflag -4 comes back with x empty.)
867 hasSolution = ~isempty(x) && all(isfinite(x));
868 if ~isempty(x) && ~hasSolution
869 warning('qrf_bas:nonFiniteSolution', ...
870 'linprog returned a non-finite solution (exitflag %d); U and the other metric fields are left unpopulated.', ...
871 exitflag);
872 end
873
874 if hasSolution
875 % Compute utilizations
876 result.U = zeros(M, 1);
877 result.e = zeros(M, max(K));
878
879 for i = 1:M
880 for ki = 1:K(i)
881 result.e(i, ki) = x(getEIdx(i, ki));
882 end
883
884 for m = 1:MR
885 for ki = 1:K(i)
886 for ni = 1:F(i)
887 idx = getP2Idx(i, ni, ki, i, ni, ki, m);
888 result.U(i) = result.U(i) + x(idx);
889 end
890 end
891 end
892 end
893 end
894
895 if verbose; fprintf('\n=== Results ===\n'); end
896 if verbose; fprintf('Objective value: %f\n', fval); end
897 if verbose; fprintf('Exit flag: %d\n', exitflag); end
898 if hasSolution && verbose
899 fprintf('\nUtilizations:\n');
900 for i = 1:M
901 fprintf(' Queue %d: U = %.6f\n', i, result.U(i));
902 end
903 fprintf('\nEffective utilizations by phase:\n');
904 for i = 1:M
905 fprintf(' Queue %d: e = [', i);
906 for ki = 1:K(i)
907 fprintf('%.6f ', result.e(i, ki));
908 end
909 fprintf(']\n');
910 end
911 end
912end
the
Definition Station.m:245
Generated by doxygen 1.9.8