1## Copyright (C) 2012, 2016, 2020 Moreno Marzolla 2## 3## This file is part of the queueing toolbox. 4## 5## The queueing toolbox is free software: you can redistribute it and/or 6## modify it under the terms of the GNU General Public License as 7## published by the Free Software Foundation, either version 3 of the 8## License, or (at your option) any later version. 9## 10## The queueing toolbox is distributed in the hope that it will be 11## useful, but WITHOUT ANY WARRANTY; without even the implied warranty 12## of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 13## General Public License for more details. 14## 15## You should have received a copy of the GNU General Public License 16## along with the queueing toolbox. If not, see <http://www.gnu.org/licenses/>. 17 18## -*- texinfo -*- 19## 20## @deftypefn {Function File} {[@var{V} @var{ch}] =} qncmvisits (@var{P}) 21## @deftypefnx {Function File} {[@var{V} @var{ch}] =} qncmvisits (@var{P}, @var{r}) 22## 23## Compute the average number of visits for the nodes of a closed multiclass network with @math{K} service centers and @math{C} customer classes. 24## 25## @strong{INPUTS} 26## 27## @table @code 28## 29## @item @var{P}(r,i,s,j) 30## probability that a 31## class @math{r} request which completed service at center @math{i} is 32## routed to center @math{j} as a class @math{s} request. Class switching 33## is allowed. 34## 35## @item @var{r}(c) 36## index of class @math{c} reference station, 37## @math{r(c) \in @{1, @dots{}, K@}}, @math{1 @leq{} c @leq{} C}. 38## The class @math{c} visit count to server @code{@var{r}(c)} 39## (@code{@var{V}(c,r(c))}) is conventionally set to 1. The reference 40## station serves two purposes: (i) its throughput is assumed to be the 41## system throughput, and (ii) a job returning to the reference station 42## is assumed to have completed one cycle. Default is to consider 43## station 1 as the reference station for all classes. 44## 45## @end table 46## 47## @strong{OUTPUTS} 48## 49## @table @code 50## 51## @item @var{V}(c,i) 52## number of visits of class @math{c} requests at center @math{i}. 53## 54## @item @var{ch}(c) 55## chain number that class @math{c} belongs 56## to. Different classes can belong to the same chain. Chains are 57## numbered sequentially starting from 1 (@math{1, 2, @dots{}}). The 58## total number of chains is @code{max(@var{ch})}. 59## 60## @end table 61## 62## @end deftypefn 63 64## Author: Moreno Marzolla <moreno.marzolla(at)unibo.it> 65## Web: http://www.moreno.marzolla.name/ 66 67function [V chains] = qncmvisits( P, r ) 68 69 if ( nargin < 1 || nargin > 2 ) 70 print_usage(); 71 endif 72 73 ndims(P) == 4 || ... 74 error("P must be a 4-dimensional matrix"); 75 76 [C, K, C2, K2] = size( P ); 77 (K == K2 && C == C2) || ... 78 error( "P must be a C*K*C*K matrix"); 79 80 if ( nargin < 2) 81 r = ones(1,C); 82 else 83 isvector(r) && length(r) == C || ... 84 error("r must be a vector with %d elements",C); 85 all( r>=1 && r<=K ) || ... 86 error("elements in r must be in the range [1, %d]",K); 87 r = r(:)'; 88 endif 89 90 ## solve the traffic equations: V(s,j) = sum_r sum_i V(r,i) * 91 ## P(r,i,s,j), for all s,j V(s,r(s)) = 1 for all s. 92 A = reshape(P,[K*C K*C])-eye(K*C); 93 b = zeros(1,K*C); 94 95 CH = __scc(reshape(P,[C*K C*K])>0); 96 nCH = max(CH); # number of chains 97 CH = reshape(CH,C,K); # CH(c,k) is the chain that class c at center k belongs to 98 99 chains = zeros(1,C); 100 101 for k=1:K 102 for c=1:C 103 if ( chains(c) == 0 ) 104 chains(c) = CH(c,k); 105 else 106 ( CH(c,k) == 0 || chains(c) == CH(c,k) ) || ... 107 error("Class %d belongs to different chains",c); 108 endif 109 endfor 110 endfor 111 112 constraints = zeros(1,nCH); # constraint(cc) = 1 iff we set a constraint for a class belonging to chain cc; we only set one constraint per chain 113 114 for c=1:C 115 cc = CH(c,r(c)); 116 if ( cc == 0 || constraints(cc) == 0 ) 117 ii = sub2ind([C K],c,r(c)); 118 A(:,ii) = 0; 119 A(ii,ii) = 1; 120 if ( cc > 0 ) ## if r(c) is not an isolated node 121 constraints(cc) = 1; 122 b(ii) = 1; 123 endif 124 endif 125 endfor 126 127 V = reshape(b/A, C, K); 128 ## Make sure that no negative values appear (sometimes, numerical 129 ## errors produce tiny negative values instead of zeros) 130 V = max(0,V); 131endfunction 132 133%!test 134%! 135%! ## Closed, multiclass network 136%! 137%! C = 2; K = 3; 138%! P = zeros(C,K,C,K); 139%! P(1,1,1,2) = 1; 140%! P(1,2,1,1) = 1; 141%! P(2,1,2,3) = 1; 142%! P(2,3,2,1) = 1; 143%! V = qncmvisits(P); 144%! for c=1:C 145%! for k=1:K 146%! assert(V(c,k), sum(sum(V .* P(:,:,c,k))), 1e-5); 147%! endfor 148%! endfor 149 150%!test 151%! 152%! ## Test multiclass network. Example from Schwetman (figure 7, page 9 of 153%! ## http://docs.lib.purdue.edu/cstech/259/ 154%! ## "Testing network-of-queues software, technical report CSD-TR 330, 155%! ## Purdue University). 156%! 157%! C = 2; K = 4; 158%! P = zeros(C,K,C,K); 159%! # class 1 routing 160%! P(1,1,1,1) = .05; 161%! P(1,1,1,2) = .45; 162%! P(1,1,1,3) = .5; 163%! P(1,2,1,1) = 1; 164%! P(1,3,1,1) = 1; 165%! # class 2 routing 166%! P(2,1,2,1) = .01; 167%! P(2,1,2,3) = .5; 168%! P(2,1,2,4) = .49; 169%! P(2,3,2,1) = 1; 170%! P(2,4,2,1) = 1; 171%! V = qncmvisits(P); 172%! for c=1:C 173%! for i=1:K 174%! assert(V(c,i), sum(sum(V .* P(:,:,c,i))), 1e-5); 175%! endfor 176%! endfor 177 178%!test 179%! 180%! ## Network with class switching. 181%! ## This is the example in figure 9 of 182%! ## Schwetman, "Implementing the Mean Value Analysis 183%! ## Algorithm fort the solution of Queueing Network Models", Technical 184%! ## Report CSD-TR-355, Department of Computer Science, Purdue Univrsity, 185%! ## Feb 15, 1982, http://docs.lib.purdue.edu/cstech/286/ 186%! 187%! C = 2; K = 3; 188%! S = [.01 .07 .10; ... 189%! .05 0.7 .10 ]; 190%! P = zeros(C,K,C,K); 191%! P(1,1,1,2) = .7; 192%! P(1,1,1,3) = .2; 193%! P(1,1,2,1) = .1; 194%! P(2,1,2,2) = .3; 195%! P(2,1,2,3) = .5; 196%! P(2,1,1,1) = .2; 197%! P(1,2,1,1) = P(1,3,1,1) = 1; 198%! P(2,2,2,1) = P(2,3,2,1) = 1; 199%! N = [3 0]; 200%! V = qncmvisits(P); 201%! VV = [10 7 2; 5 1.5 2.5]; # result given in Schwetman; our function computes something different, but that's ok since visit counts are actually ratios 202%! assert( V ./ repmat(V(:,1),1,K), VV ./ repmat(VV(:,1),1,K), 1e-5 ); 203 204%!test 205%! 206%! ## two disjoint classes: must produce two disjoing chains 207%! 208%! C = 2; K = 3; 209%! P = zeros(C,K,C,K); 210%! P(1,1,1,2) = 1; 211%! P(1,2,1,1) = 1; 212%! P(2,1,2,3) = 1; 213%! P(2,3,2,1) = 1; 214%! [nc r] = qncmvisits(P); 215%! assert( r(1) != r(2) ); 216 217%!test 218%! 219%! ## two classes, one chain 220%! 221%! C = 2; K = 3; 222%! P = zeros(C,K,C,K); 223%! P(1,1,1,2) = .5; 224%! P(1,2,2,1) = 1; 225%! P(2,1,2,3) = .5; 226%! P(2,3,1,1) = 1; 227%! [nc r] = qncmvisits(P); 228%! assert( r(1) == r(2) ); 229 230%!test 231%! 232%! ## a "Moebius strip". Note that this configuration is invalid, and 233%! ## therefore our algorithm must raise an error. This is because this 234%! ## network has two chains, but both chains contain both classes 235%! 236%! C = 2; K = 2; 237%! P = zeros(C,K,C,K); 238%! P(1,1,2,2) = 1; 239%! P(2,2,1,1) = 1; 240%! P(2,1,1,2) = 1; 241%! P(1,2,2,1) = 1; 242%! fail( "qncmvisits(P)", "different"); 243 244%!test 245%! 246%! ## Network with two classes representing independent chains. 247%! ## This is example in figure 8 of 248%! ## Schwetman, "Implementing the Mean Value Analysis 249%! ## Algorithm fort the solution of Queueing Network Models", Technical 250%! ## Report CSD-TR-355, Department of Computer Science, Purdue Univrsity, 251%! ## Feb 15, 1982, http://docs.lib.purdue.edu/cstech/286/ 252%! 253%! C = 2; K = 2; 254%! P = zeros(C,K,C,K); 255%! P(1,1,1,3) = P(1,3,1,1) = 1; 256%! P(2,2,2,3) = P(2,3,2,2) = 1; 257%! V = qncmvisits(P,[1,2]); 258%! assert( V, [1 0 1; 0 1 1], 1e-5 ); 259 260%!test 261%! C = 2; 262%! K = 3; 263%! P = zeros(C,K,C,K); 264%! P(1,1,1,2) = 1; 265%! P(1,2,1,3) = 1; 266%! P(1,3,2,2) = 1; 267%! P(2,2,1,1) = 1; 268%! [V ch] = qncmvisits(P); 269%! assert( ch, [1 1] ); 270 271## The following transition probability matrix is not well formed: note 272## that there is an outgoing transition from center 1, class 1 but not 273## incoming transition. 274%!test 275%! C = 2; 276%! K = 3; 277%! P = zeros(C,K,C,K); 278%! P(1,1,1,2) = 1; 279%! P(1,2,1,3) = 1; 280%! P(1,3,2,2) = 1; 281%! P(2,2,2,1) = 1; 282%! P(2,1,1,2) = 1; 283%! [V ch] = qncmvisits(P); 284%! assert( ch, [1 1] ); 285 286## compute strongly connected components using Kosaraju's algorithm, 287## which requires two DFS visits. A better solution would be to use 288## Tarjan's algorithm. 289## 290## In this implementation, an isolated node without self loops will NOT 291## belong to any SCC. Although this is not formally correct from the 292## graph theoretic point of view, it is necessary to compute chains 293## correctly. 294function s = __scc(G) 295 assert(issquare(G)); 296 N = rows(G); 297 GF = (G>0); 298 GB = (G'>0); 299 s = zeros(N,1); 300 c=1; 301 for n=1:N 302 if (s(n) == 0) 303 fw = __dfs(GF,n); 304 bw = __dfs(GB,n); 305 r = (fw & bw); 306 if (any(r)) 307 s(r) = c++; 308 endif 309 endif 310 endfor 311endfunction 312 313## Executes a dfs visit on graph G, starting from source node s 314function v = __dfs(G, s) 315 assert( issquare(G) ); 316 N = rows(G); 317 v = stack = zeros(1,N); ## v(i) == 1 iff node i has been visited 318 q = 1; # first empty slot in queue 319 stack(q++) = s; 320 ## Note: node s is NOT marked as visited; it will be marked as visited 321 ## only if we visit it again. This is necessary to ensure that 322 ## isolated nodes without self loops will not belong to any SCC. 323 while( q>1 ) 324 n = stack(--q); 325 ## explore neighbors of n: all f in G(n,:) such that v(f) == 0 326 327 ## The following instruction is equivalent to: 328 ## for f=find(G(n,:)) 329 ## if ( v(f) == 0 ) 330 for f = find ( G(n,:) & (v==0) ) 331 stack(q++) = f; 332 v(f) = 1; 333 endfor 334 endwhile 335endfunction 336 337