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% objective - (optional) 'U1min' (default), 'U1max', or queue index 1..M
29% sense - (optional) 'min' (default) or 'max'
30%
31% Outputs:
32% result - Structure with results:
33% .U - [M x 1] Utilization of each queue
34% .e - [M x max(K)] Effective utilization by phase
35% x - Raw solution vector
36% fval - Objective function value
37% exitflag - Solver exit flag
38
39 if nargin < 2 || isempty(objective)
40 objective = 'U1min';
41 end
42 if nargin < 3 || isempty(sense)
43 sense = 'min';
44 end
45
46 % Extract parameters
47 M = params.M;
48 N = params.N;
49 f = params.f; % finite capacity queue index
50 F = params.F(:);
51 K = params.K(:);
52 mu = params.mu;
53 v = params.v;
54 r = params.r;
55 MR = params.MR;
56 BB = params.BB;
57 MM = params.MM;
58 ZZ = params.ZZ(:);
59 ZM = params.ZM;
60 MM1 = params.MM1;
61 if isfield(params, 'verbose')
62 verbose = params.verbose;
63 else
64 verbose = true;
65 end
66
67 % Compute transition rates q(i,j,k,h)
68 % q{i,j} is a K(i) x K(i) array (not load-dependent for BAS)
69 q = cell(M, M);
70 for i = 1:M
71 for j = 1:M
72 q{i,j} = zeros(K(i), K(i));
73 for ki = 1:K(i)
74 for hi = 1:K(i)
75 if j ~= i
76 q{i,j}(ki, hi) = r(i,j) * mu{i}(ki, hi);
77 else
78 q{i,j}(ki, hi) = v{i}(ki, hi) + r(i,i) * mu{i}(ki, hi);
79 end
80 end
81 end
82 end
83 end
84
85 %% Build variable indexing
86 % p2(j, nj, kj, i, ni, hi, m) for j in 1:M, nj in 0:N, kj in 1:K(j),
87 % i in 1:M, ni in 0:N, hi in 1:K(i), m in 1:MR
88 % e(i, ki) for i in 1:M, ki in 1:K(i)
89
90 if verbose; fprintf('Building variable index map...\n'); end
91 varCount = 0;
92
93 % p2 variables
94 p2idx = cell(M, 1);
95 for j = 1:M
96 p2idx{j} = cell(N+1, K(j), M, MR);
97 for nj = 0:N
98 for kj = 1:K(j)
99 for i = 1:M
100 for m = 1:MR
101 p2idx{j}{nj+1, kj, i, m} = zeros(N+1, K(i));
102 for ni = 0:N
103 for hi = 1:K(i)
104 varCount = varCount + 1;
105 p2idx{j}{nj+1, kj, i, m}(ni+1, hi) = varCount;
106 end
107 end
108 end
109 end
110 end
111 end
112 end
113
114 % e variables
115 eidx = cell(M, 1);
116 for i = 1:M
117 eidx{i} = zeros(K(i), 1);
118 for ki = 1:K(i)
119 varCount = varCount + 1;
120 eidx{i}(ki) = varCount;
121 end
122 end
123
124 nVars = varCount;
125 if verbose; fprintf('Total variables: %d\n', nVars); end
126
127 %% Build constraints
128 Aeq = [];
129 beq = [];
130 Aineq = [];
131 bineq = [];
132
133 % Helper functions
134 getP2Idx = @(j, nj, kj, i, ni, hi, m) p2idx{j}{nj+1, kj, i, m}(ni+1, hi);
135 getEIdx = @(i, ki) eidx{i}(ki);
136
137 if verbose; fprintf('Building constraints...\n'); end
138
139 %% Initialize bounds
140 lb = zeros(nVars, 1);
141 ub = inf(nVars, 1);
142
143 %% ZERO constraints - fix infeasible states
144 if verbose; fprintf(' ZERO constraints...\n'); end
145 for j = 1:M
146 for nj = 0:N
147 for kj = 1:K(j)
148 for i = 1:M
149 for ni = 0:N
150 for hi = 1:K(i)
151 for m = 1:MR
152 idx = getP2Idx(j, nj, kj, i, ni, hi, m);
153
154 % ZERO1: i==j, nj==ni, h<>k
155 if i == j && nj == ni && hi ~= kj
156 ub(idx) = 0;
157 end
158
159 % ZERO2: i==j, nj<>ni
160 if i == j && nj ~= ni
161 ub(idx) = 0;
162 end
163
164 % ZERO3: i<>j, nj+ni > N
165 if i ~= j && nj + ni > N
166 ub(idx) = 0;
167 end
168
169 % ZERO6: nj > F(j)
170 if nj > F(j)
171 ub(idx) = 0;
172 end
173
174 % ZERO5: BB(m,j)==1 and nj==0
175 if m >= 2 && BB(m, j) == 1 && nj == 0
176 ub(idx) = 0;
177 end
178
179 % ZERO7: BB(m,j)==1 and i<>j and i<>f and ni+nj+F(f)>N
180 if m >= 2 && BB(m, j) == 1 && i ~= j && i ~= f && ni + nj + F(f) > N
181 ub(idx) = 0;
182 end
183
184 % ZERO8: finite queue not at capacity in blocking config
185 if j == f && nj >= 1 && nj <= F(f)-1 && m >= 2
186 ub(idx) = 0;
187 end
188 end
189 end
190 end
191 end
192 end
193 end
194 end
195
196 % ZERO4: For m>=2 and j<>f, p2(j,nj,k,f,nf,h,m)=0 when nf < F(f)
197 for j = 1:M
198 if j == f
199 continue;
200 end
201 for nj = 0:N
202 for kj = 1:K(j)
203 for m = 2:MR
204 for nf = 0:(F(f)-1)
205 for hf = 1:K(f)
206 idx = getP2Idx(j, nj, kj, f, nf, hf, m);
207 ub(idx) = 0;
208 end
209 end
210 end
211 end
212 end
213 end
214
215 %% ONE: Normalization
216 if verbose; fprintf(' ONE constraints...\n'); end
217 for j = 1:M
218 row = zeros(1, nVars);
219 for nj = 0:N
220 for kj = 1:K(j)
221 for m = 1:MR
222 idx = getP2Idx(j, nj, kj, j, nj, kj, m);
223 row(idx) = 1;
224 end
225 end
226 end
227 Aeq = [Aeq; row];
228 beq = [beq; 1];
229 end
230
231 %% SYMMETRY
232 if verbose; fprintf(' SYMMETRY constraints...\n'); end
233 for j = 1:M
234 for nj = 0:min(N, F(j))
235 for kj = 1:K(j)
236 for i = 1:M
237 if i <= j
238 continue;
239 end
240 for ni = 0:min(N, F(i))
241 if i ~= j && nj + ni > N
242 continue;
243 end
244 for hi = 1:K(i)
245 for m = 1:MR
246 idx1 = getP2Idx(j, nj, kj, i, ni, hi, m);
247 idx2 = getP2Idx(i, ni, hi, j, nj, kj, m);
248 if ub(idx1) == 0 && ub(idx2) == 0
249 continue;
250 end
251 if idx1 ~= idx2
252 row = zeros(1, nVars);
253 row(idx1) = 1;
254 row(idx2) = -1;
255 Aeq = [Aeq; row];
256 beq = [beq; 0];
257 end
258 end
259 end
260 end
261 end
262 end
263 end
264 end
265
266 %% MARGINALS
267 if verbose; fprintf(' MARGINALS constraints...\n'); end
268 for j = 1:M
269 for kj = 1:K(j)
270 for nj = 0:min(N, F(j))
271 for i = 1:M
272 if i == j
273 continue;
274 end
275 for m = 1:MR
276 row = zeros(1, nVars);
277 idx_diag = getP2Idx(j, nj, kj, j, nj, kj, m);
278 row(idx_diag) = 1;
279 for ni = 0:min(N-nj, F(i))
280 for hi = 1:K(i)
281 idx = getP2Idx(j, nj, kj, i, ni, hi, m);
282 row(idx) = row(idx) - 1;
283 end
284 end
285 Aeq = [Aeq; row];
286 beq = [beq; 0];
287 end
288 end
289 end
290 end
291 end
292
293 %% UEFF: e(i,ki) = sum of p2 where queue i is not blocked
294 if verbose; fprintf(' UEFF constraints...\n'); end
295 for i = 1:M
296 for ki = 1:K(i)
297 row = zeros(1, nVars);
298 idx_e = getEIdx(i, ki);
299 row(idx_e) = -1;
300 for j = 1:M
301 for nj = 0:min(N, F(j))
302 for kj = 1:K(j)
303 for m = 1:MR
304 if BB(m, i) == 0
305 for ni = 1:min(N, F(i))
306 idx = getP2Idx(j, nj, kj, i, ni, ki, m);
307 row(idx) = row(idx) + 1;
308 end
309 end
310 end
311 end
312 end
313 end
314 Aeq = [Aeq; row];
315 beq = [beq; 0];
316 end
317 end
318
319 %% THM1: Phase balance (Theorem 1)
320 % 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)
321 if verbose; fprintf(' THM1 (Phase balance) constraints...\n'); end
322 for i = 1:M
323 for ki = 1:K(i)
324 row = zeros(1, nVars);
325 % LHS
326 for j = 1:M
327 for hi = 1:K(i)
328 if j ~= i || hi ~= ki
329 idx_e = getEIdx(i, ki);
330 row(idx_e) = row(idx_e) + q{i,j}(ki, hi);
331 end
332 end
333 end
334 % RHS (subtract)
335 for j = 1:M
336 for hi = 1:K(i)
337 if j ~= i || hi ~= ki
338 idx_e = getEIdx(i, hi);
339 row(idx_e) = row(idx_e) - q{i,j}(hi, ki);
340 end
341 end
342 end
343 Aeq = [Aeq; row];
344 beq = [beq; 0];
345 end
346 end
347
348 %% THM2: Population constraint (Theorem 2)
349 if verbose; fprintf(' THM2 (Population) constraints...\n'); end
350 for j = 1:M
351 for kj = 1:K(j)
352 for nj = 0:F(j)
353 for m = 1:MR
354 row = zeros(1, nVars);
355 % RHS: -N * p2(j,nj,kj,j,nj,kj,m)
356 idx_diag = getP2Idx(j, nj, kj, j, nj, kj, m);
357 row(idx_diag) = -N;
358 % LHS: sum
359 for i = 1:M
360 for ni = 1:F(i)
361 for ki = 1:K(i)
362 idx = getP2Idx(j, nj, kj, i, ni, ki, m);
363 row(idx) = row(idx) + ni;
364 end
365 end
366 end
367 Aeq = [Aeq; row];
368 beq = [beq; 0];
369 end
370 end
371 end
372 end
373
374 %% COR1: Second moment constraint (Corollary to Theorem 2)
375 % sum_{m,i,j,nj,ni,ki,kj} ni*nj*p2(j,nj,kj,i,ni,ki,m) = N^2
376 if verbose; fprintf(' COR1 (Second moment) constraint...\n'); end
377 row = zeros(1, nVars);
378 for m = 1:MR
379 for i = 1:M
380 for j = 1:M
381 for nj = 1:F(j)
382 for ni = 1:F(i)
383 for ki = 1:K(i)
384 for kj = 1:K(j)
385 idx = getP2Idx(j, nj, kj, i, ni, ki, m);
386 row(idx) = row(idx) + ni * nj;
387 end
388 end
389 end
390 end
391 end
392 end
393 end
394 Aeq = [Aeq; row];
395 beq = [beq; N^2];
396
397 %% THM30: Marginal balance for ni=0 (per phase), i<>f
398 if verbose; fprintf(' THM30 (Marginal balance ni=0) constraints...\n'); end
399 for i = 1:M
400 if i == f
401 continue;
402 end
403 for ui = 1:K(i)
404 row = zeros(1, nVars);
405 % LHS: arrivals from j<>i,j<>f with BB(m,j)==0
406 for j = 1:M
407 if j == i || j == f
408 continue;
409 end
410 for nj = 1:F(j)
411 for kj = 1:K(j)
412 for hj = 1:K(j)
413 for m = 1:MR
414 if BB(m, j) == 0
415 idx = getP2Idx(j, nj, kj, i, 0, ui, m);
416 row(idx) = row(idx) + q{j,i}(kj, hj);
417 end
418 end
419 end
420 end
421 end
422 end
423 % LHS: arrivals from j==f with MM(m,1)<>i
424 for nj = 1:F(f)
425 for kj = 1:K(f)
426 for hj = 1:K(f)
427 for m = 1:MR
428 if MM(m, 1) ~= i
429 idx = getP2Idx(f, nj, kj, i, 0, ui, m);
430 row(idx) = row(idx) + q{f,i}(kj, hj);
431 end
432 end
433 end
434 end
435 end
436
437 % RHS: departures from i at ni=1 to j<>i,j<>f with BB(m,i)==0
438 for j = 1:M
439 if j == i || j == f
440 continue;
441 end
442 for nj = 0:F(j)
443 for ki = 1:K(i)
444 for hj = 1:K(j)
445 for m = 1:MR
446 if BB(m, i) == 0
447 idx = getP2Idx(j, nj, hj, i, 1, ki, m);
448 row(idx) = row(idx) - q{i,j}(ki, ui);
449 end
450 end
451 end
452 end
453 end
454 end
455 % RHS: departures to j==f with BB(m,i)==0 and nj<F(f)
456 for nj = 0:(F(f)-1)
457 for ki = 1:K(i)
458 for hj = 1:K(f)
459 for m = 1:MR
460 if BB(m, i) == 0
461 idx = getP2Idx(f, nj, hj, i, 1, ki, m);
462 row(idx) = row(idx) - q{i,f}(ki, ui);
463 end
464 end
465 end
466 end
467 end
468 % RHS: unblocking when BB(m,i)==1 and MM(m,1)==i
469 for m = 1:MR
470 if BB(m, i) == 1 && MM(m, 1) == i
471 for kf = 1:K(f)
472 for pf = 1:K(f)
473 for w = 1:M
474 if w ~= f && w ~= i
475 idx = getP2Idx(f, F(f), kf, i, 1, ui, m);
476 row(idx) = row(idx) - q{f,w}(kf, pf);
477 end
478 end
479 end
480 end
481 end
482 end
483
484 Aeq = [Aeq; row];
485 beq = [beq; 0];
486 end
487 end
488
489 %% THM3: Marginal balance for ni in 1:F(i)-1, i<>f
490 if verbose; fprintf(' THM3 (Marginal balance) constraints...\n'); end
491 for i = 1:M
492 if i == f
493 continue;
494 end
495 for ni = 1:(F(i)-1)
496 row = zeros(1, nVars);
497 % LHS: arrivals from j<>i,j<>f with BB(m,j)==0
498 for j = 1:M
499 if j == i || j == f
500 continue;
501 end
502 for nj = 1:F(j)
503 for kj = 1:K(j)
504 for hj = 1:K(j)
505 for ui = 1:K(i)
506 for m = 1:MR
507 if BB(m, j) == 0
508 idx = getP2Idx(j, nj, kj, i, ni, ui, m);
509 row(idx) = row(idx) + q{j,i}(kj, hj);
510 end
511 end
512 end
513 end
514 end
515 end
516 end
517 % LHS: arrivals from j==f with MM(m,1)<>i
518 for nj = 1:F(f)
519 for kj = 1:K(f)
520 for hj = 1:K(f)
521 for ui = 1:K(i)
522 for m = 1:MR
523 if MM(m, 1) ~= i
524 idx = getP2Idx(f, nj, kj, i, ni, ui, m);
525 row(idx) = row(idx) + q{f,i}(kj, hj);
526 end
527 end
528 end
529 end
530 end
531 end
532
533 % RHS: departures from i at ni+1 to j<>i,j<>f with BB(m,i)==0
534 for j = 1:M
535 if j == i || j == f
536 continue;
537 end
538 for nj = 0:F(j)
539 for ki = 1:K(i)
540 for hi = 1:K(i)
541 for uj = 1:K(j)
542 for m = 1:MR
543 if BB(m, i) == 0
544 idx = getP2Idx(j, nj, uj, i, ni+1, ki, m);
545 row(idx) = row(idx) - q{i,j}(ki, hi);
546 end
547 end
548 end
549 end
550 end
551 end
552 end
553 % RHS: departures to j==f with BB(m,i)==0 and nj<F(f)
554 for nj = 0:(F(f)-1)
555 for ki = 1:K(i)
556 for hi = 1:K(i)
557 for uj = 1:K(f)
558 for m = 1:MR
559 if BB(m, i) == 0
560 idx = getP2Idx(f, nj, uj, i, ni+1, ki, m);
561 row(idx) = row(idx) - q{i,f}(ki, hi);
562 end
563 end
564 end
565 end
566 end
567 end
568 % RHS: unblocking when BB(m,i)==1 and MM(m,1)==i
569 for m = 1:MR
570 if BB(m, i) == 1 && MM(m, 1) == i
571 for ki = 1:K(i)
572 for kf = 1:K(f)
573 for pf = 1:K(f)
574 for w = 1:M
575 if w ~= f && w ~= i
576 idx = getP2Idx(f, F(f), kf, i, ni+1, ki, m);
577 row(idx) = row(idx) - q{f,w}(kf, pf);
578 end
579 end
580 end
581 end
582 end
583 end
584 end
585
586 Aeq = [Aeq; row];
587 beq = [beq; 0];
588 end
589 end
590
591 %% THM3f: Marginal balance for i==f
592 if verbose; fprintf(' THM3f (Marginal balance for finite queue) constraints...\n'); end
593 for ni = 0:(F(f)-1)
594 row = zeros(1, nVars);
595 % LHS: arrivals from j<>f with BB(m,j)==0 and ni<F(f)
596 if ni < F(f)
597 for j = 1:M
598 if j == f
599 continue;
600 end
601 for nj = 1:F(j)
602 for kj = 1:K(j)
603 for hj = 1:K(j)
604 for uf = 1:K(f)
605 for m = 1:MR
606 if BB(m, j) == 0
607 idx = getP2Idx(j, nj, kj, f, ni, uf, m);
608 row(idx) = row(idx) + q{j,f}(kj, hj);
609 end
610 end
611 end
612 end
613 end
614 end
615 end
616 end
617
618 % RHS: departures from f at ni+1 to j<>f (only m=1, no blocking)
619 for j = 1:M
620 if j == f
621 continue;
622 end
623 for nj = 0:F(j)
624 for kf = 1:K(f)
625 for hf = 1:K(f)
626 for uj = 1:K(j)
627 if ni < F(f)
628 idx = getP2Idx(j, nj, uj, f, ni+1, kf, 1);
629 row(idx) = row(idx) - q{f,j}(kf, hf);
630 end
631 end
632 end
633 end
634 end
635 end
636
637 Aeq = [Aeq; row];
638 beq = [beq; 0];
639 end
640
641 %% THM3I: Blocking depth balance (Theorem 4)
642 if verbose; fprintf(' THM3I (Blocking depth balance) constraints...\n'); end
643 for z = 0:(ZM-1)
644 row = zeros(1, nVars);
645 % LHS: arrivals to f at F(f) from j<>f with BB(m,j)==0 and ZZ(m)==z
646 for j = 1:M
647 if j == f
648 continue;
649 end
650 for nj = 1:F(j)
651 for kj = 1:K(j)
652 for hj = 1:K(j)
653 for uf = 1:K(f)
654 for m = 1:MR
655 if BB(m, j) == 0 && ZZ(m) == z
656 idx = getP2Idx(j, nj, kj, f, F(f), uf, m);
657 row(idx) = row(idx) + q{j,f}(kj, hj);
658 end
659 end
660 end
661 end
662 end
663 end
664 end
665 % RHS: departures from f to j<>f with ZZ(m)==z+1
666 for j = 1:M
667 if j == f
668 continue;
669 end
670 for nj = 0:F(j)
671 for kf = 1:K(f)
672 for hf = 1:K(f)
673 for uj = 1:K(j)
674 for m = 1:MR
675 if ZZ(m) == z+1
676 idx = getP2Idx(j, nj, uj, f, F(f), kf, m);
677 row(idx) = row(idx) - q{f,j}(kf, hf);
678 end
679 end
680 end
681 end
682 end
683 end
684 end
685
686 Aeq = [Aeq; row];
687 beq = [beq; 0];
688 end
689
690 %% THM3L: Maximum blocking depth constraint
691 if verbose; fprintf(' THM3L (Max blocking depth) constraints...\n'); end
692 for m = 1:MR
693 if ZZ(m) ~= ZM - 1
694 continue;
695 end
696 row = zeros(1, nVars);
697 % LHS: arrivals from j<>f with BB(m,j)==0 and MM1(m,j)>0
698 for j = 1:M
699 if j == f || BB(m, j) ~= 0 || MM1(m, j) <= 0
700 continue;
701 end
702 for nj = 1:F(j)
703 for kj = 1:K(j)
704 for hj = 1:K(j)
705 for uf = 1:K(f)
706 idx = getP2Idx(j, nj, kj, f, F(f), uf, m);
707 row(idx) = row(idx) + q{j,f}(kj, hj);
708 end
709 end
710 end
711 end
712 end
713 % RHS: uses MM1(m,j) to index into blocking configuration
714 for j = 1:M
715 if j == f || BB(m, j) ~= 0 || MM1(m, j) <= 0
716 continue;
717 end
718 mp = MM1(m, j); % blocking configuration index
719 for kf = 1:K(f)
720 for uf = 1:K(f)
721 for w = 1:M
722 if w ~= f
723 idx = getP2Idx(f, F(f), kf, f, F(f), kf, mp);
724 row(idx) = row(idx) - q{f,w}(kf, uf);
725 end
726 end
727 end
728 end
729 end
730
731 Aeq = [Aeq; row];
732 beq = [beq; 0];
733 end
734
735 %% THM4: Queue-length bound inequality (Theorem 5)
736 if verbose; fprintf(' THM4 (Queue-length bound) constraints...\n'); end
737 for j = 1:M
738 for kj = 1:K(j)
739 for i = 1:M
740 for m = 1:MR
741 row = zeros(1, nVars);
742 % LHS: sum_t sum_ht sum_nj sum_nt nt * p2
743 for t = 1:M
744 for ht = 1:K(t)
745 for nj = 0:F(j)
746 for nt = 1:F(t)
747 idx = getP2Idx(j, nj, kj, t, nt, ht, m);
748 row(idx) = row(idx) + nt;
749 end
750 end
751 end
752 end
753 % RHS: -N * sum
754 for hi = 1:K(i)
755 for nj = 0:F(j)
756 for ni = 1:F(i)
757 idx = getP2Idx(j, nj, kj, i, ni, hi, m);
758 row(idx) = row(idx) - N;
759 end
760 end
761 end
762 Aineq = [Aineq; -row];
763 bineq = [bineq; 0];
764 end
765 end
766 end
767 end
768
769 %% Build objective function
770 if verbose; fprintf('Building objective function...\n'); end
771 c = zeros(nVars, 1);
772
773 if ischar(objective)
774 if strcmp(objective, 'U1min') || strcmp(objective, 'U1max')
775 targetQueue = 1;
776 else
777 error('Unknown objective: %s', objective);
778 end
779 else
780 targetQueue = objective;
781 end
782
783 % Utilization = sum over m, k, n of p2(i,n,k,i,n,k,m)
784 for m = 1:MR
785 for ki = 1:K(targetQueue)
786 for ni = 1:F(targetQueue)
787 idx = getP2Idx(targetQueue, ni, ki, targetQueue, ni, ki, m);
788 c(idx) = 1;
789 end
790 end
791 end
792
793 if strcmp(sense, 'max')
794 c = -c;
795 end
796
797 %% Solve LP
798 if verbose; fprintf('Solving LP with %d variables and %d equality + %d inequality constraints...\n', ...
799 nVars, size(Aeq, 1), size(Aineq, 1)); end
800
801 if verbose
802 options = optimoptions('linprog', 'Display', 'final', 'Algorithm', 'interior-point');
803 else
804 options = optimoptions('linprog', 'Display', 'off', 'Algorithm', 'interior-point');
805 end
806
807 [x, fval, exitflag] = linprog(c, Aineq, bineq, Aeq, beq, lb, ub, options);
808
809 if strcmp(sense, 'max')
810 fval = -fval;
811 end
812
813 %% Extract results
814 result = struct();
815 result.objective = fval;
816 result.exitflag = exitflag;
817
818 if exitflag > 0
819 % Compute utilizations
820 result.U = zeros(M, 1);
821 result.e = zeros(M, max(K));
822
823 for i = 1:M
824 for ki = 1:K(i)
825 result.e(i, ki) = x(getEIdx(i, ki));
826 end
827
828 for m = 1:MR
829 for ki = 1:K(i)
830 for ni = 1:F(i)
831 idx = getP2Idx(i, ni, ki, i, ni, ki, m);
832 result.U(i) = result.U(i) + x(idx);
833 end
834 end
835 end
836 end
837 end
838
839 if verbose; fprintf('\n=== Results ===\n'); end
840 if verbose; fprintf('Objective value: %f\n', fval); end
841 if verbose; fprintf('Exit flag: %d\n', exitflag); end
842 if exitflag > 0 && verbose
843 fprintf('\nUtilizations:\n');
844 for i = 1:M
845 fprintf(' Queue %d: U = %.6f\n', i, result.U(i));
846 end
847 fprintf('\nEffective utilizations by phase:\n');
848 for i = 1:M
849 fprintf(' Queue %d: e = [', i);
850 for ki = 1:K(i)
851 fprintf('%.6f ', result.e(i, ki));
852 end
853 fprintf(']\n');
854 end
855 end
856end
Generated by doxygen 1.9.8