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