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} =} qncsvisits (@var{P}) 21## @deftypefnx {Function File} {@var{V} =} qncsvisits (@var{P}, @var{r}) 22## 23## Compute the mean number of visits to the service centers of a 24## single class, closed network with @math{K} service centers. 25## 26## @strong{INPUTS} 27## 28## @table @code 29## 30## @item @var{P}(i,j) 31## probability that a request which completed service at center 32## @math{i} is routed to center @math{j} (@math{K \times K} matrix). 33## For closed networks it must hold that @code{sum(@var{P},2)==1}. The 34## routing graph must be strongly connected, meaning that each node 35## must be reachable from every other node. 36## 37## @item @var{r} 38## Index of the reference station, @math{r \in @{1, @dots{}, K@}}; 39## Default @code{@var{r}=1}. The traffic equations are solved by 40## imposing the condition @code{@var{V}(r) = 1}. A request returning to 41## the reference station completes its activity cycle. 42## 43## @end table 44## 45## @strong{OUTPUTS} 46## 47## @table @code 48## 49## @item @var{V}(k) 50## average number of visits to service center @math{k}, assuming 51## @math{r} as the reference station. 52## 53## @end table 54## 55## @end deftypefn 56 57## Author: Moreno Marzolla <moreno.marzolla(at)unibo.it> 58## Web: http://www.moreno.marzolla.name/ 59 60function V = qncsvisits( P, r ) 61 62 if ( nargin < 1 || nargin > 2 ) 63 print_usage(); 64 endif 65 66 issquare(P) || ... 67 error("P must be a square matrix"); 68 69 [res err] = dtmcchkP(P); 70 (res>0) || ... 71 error( "invalid transition probability matrix P" ); 72 73 K = rows(P); 74 75 if ( nargin < 2 ) 76 r = 1; 77 else 78 isscalar(r) || ... 79 error("r must be a scalar"); 80 81 (r>=1 && r<=K) || ... 82 error("r must be an integer in the range [1, %d]",K); 83 endif 84 85 V = zeros(size(P)); 86 A = P-eye(K); 87 b = zeros(1,K); 88 A(:,r) = 0; A(r,r) = 1; 89 b(r) = 1; 90 V = b/A; 91 ## Make sure that no negative values appear (sometimes, numerical 92 ## errors produce tiny negative values instead of zeros) 93 V = max(0,V); 94endfunction 95%!test 96%! P = [-1 0; 0 0]; 97%! fail( "qncsvisits(P)", "invalid" ); 98%! P = [1 0; 0.5 0]; 99%! fail( "qncsvisits(P)", "invalid" ); 100%! P = [1 2 3; 1 2 3]; 101%! fail( "qncsvisits(P)", "square" ); 102%! P = [0 1; 1 0]; 103%! fail( "qncsvisits(P,0)", "range" ); 104%! fail( "qncsvisits(P,3)", "range" ); 105 106%!test 107%! 108%! ## Closed, single class network 109%! 110%! P = [0 0.3 0.7; 1 0 0; 1 0 0]; 111%! V = qncsvisits(P); 112%! assert( V*P,V,1e-5 ); 113%! assert( V, [1 0.3 0.7], 1e-5 ); 114 115%!test 116%! 117%! ## Test tolerance of the qncsvisits() function. 118%! ## This test builds transition probability matrices and tries 119%! ## to compute the visit counts on them. 120%! 121%! for k=[5, 10, 20, 50] 122%! P = reshape(1:k^2, k, k); 123%! P = P ./ repmat(sum(P,2),1,k); 124%! V = qncsvisits(P); 125%! assert( V*P, V, 1e-5 ); 126%! endfor 127 128%!demo 129%! P = [0 0.3 0.7; ... 130%! 1 0 0 ; ... 131%! 1 0 0 ]; 132%! V = qncsvisits(P) 133 134