1## Copyright (C) 2009, 2010, 2011, 2012, 2016, 2018 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{U}, @var{R}, @var{Q}, @var{X}] =} qncsmvald (@var{N}, @var{S}, @var{V}) 21## @deftypefnx {Function File} {[@var{U}, @var{R}, @var{Q}, @var{X}] =} qncsmvald (@var{N}, @var{S}, @var{V}, @var{Z}) 22## 23## @cindex Mean Value Analysys (MVA) 24## @cindex MVA 25## @cindex closed network, single class 26## @cindex load-dependent service center 27## 28## Mean Value Analysis algorithm for closed, single class queueing 29## networks with @math{K} service centers and load-dependent service 30## times. This function supports FCFS, LCFS-PR, PS and IS nodes. For 31## networks with only fixed-rate centers and multiple-server 32## nodes, the function @code{qncsmva} is more efficient. 33## 34## @strong{INPUTS} 35## 36## @table @code 37## 38## @item @var{N} 39## Population size (number of requests in the system, @code{@var{N} @geq{} 0}). 40## If @code{@var{N} == 0}, this function returns @code{@var{U} = @var{R} = @var{Q} = @var{X} = 0} 41## 42## @item @var{S}(k,n) 43## mean service time at center @math{k} 44## where there are @math{n} requests, @math{1 @leq{} n 45## @leq{} N}. @code{@var{S}(k,n)} @math{= 1 / \mu_{k}(n)}, 46## where @math{\mu_{k}(n)} is the service rate of center @math{k} 47## when there are @math{n} requests. 48## 49## @item @var{V}(k) 50## average number of visits to service center @math{k} (@code{@var{V}(k) @geq{} 0}). 51## 52## @item @var{Z} 53## external delay ("think time", @code{@var{Z} @geq{} 0}); default 0. 54## 55## @end table 56## 57## @strong{OUTPUTS} 58## 59## @table @code 60## 61## @item @var{U}(k) 62## utilization of service center @math{k}. The 63## utilization is defined as the probability that service center 64## @math{k} is not empty, that is, @math{U_k = 1-\pi_k(0)} where 65## @math{\pi_k(0)} is the steady-state probability that there are 0 66## jobs at service center @math{k}. 67## 68## @item @var{R}(k) 69## response time on service center @math{k}. 70## 71## @item @var{Q}(k) 72## average number of requests in service center @math{k}. 73## 74## @item @var{X}(k) 75## throughput of service center @math{k}. 76## 77## @end table 78## 79## @strong{NOTES} 80## 81## In presence of load-dependent servers, the MVA algorithm is known 82## to be numerically unstable. Generally this problem manifests itself 83## as negative response times or utilization. 84## 85## @strong{REFERENCES} 86## 87## @itemize 88## @item 89## M. Reiser and S. S. Lavenberg, @cite{Mean-Value Analysis of Closed 90## Multichain Queuing Networks}, Journal of the ACM, vol. 27, n. 2, 91## April 1980, pp. 313--322. @uref{http://doi.acm.org/10.1145/322186.322195, 10.1145/322186.322195} 92## @end itemize 93## 94## This implementation is described in G. Bolch, S. Greiner, H. de Meer 95## and K. Trivedi, @cite{Queueing Networks and Markov Chains: Modeling 96## and Performance Evaluation with Computer Science Applications}, Wiley, 97## 1998, Section 8.2.4.1, ``Networks with Load-Dependent Service: Closed 98## Networks''. 99## 100## @seealso{qncsmva} 101## 102## @end deftypefn 103 104## Author: Moreno Marzolla <moreno.marzolla(at)unibo.it> 105## Web: http://www.moreno.marzolla.name/ 106 107function [U R Q X] = qncsmvald( N, S, V, Z ) 108 109 if ( nargin < 3 || nargin > 4 ) 110 print_usage(); 111 endif 112 113 isvector(V) && all(V>=0) || ... 114 error( "V must be a vector >= 0" ); 115 V = V(:)'; # make V a row vector 116 K = length(V); # Number of servers 117 isscalar(N) && N >= 0 || ... 118 error( "N must be >= 0" ); 119 ( ismatrix(S) && rows(S) == K && columns(S) >= N ) || ... 120 error( "S size mismatch: is %dx%d, should be %dx%d", rows(S), columns(S), K, N ); 121 all(S(:)>=0) || ... 122 error( "S must be >= 0" ); 123 124 if ( nargin < 4 ) 125 Z = 0; 126 else 127 isscalar(Z) && Z>=0 || ... 128 error( "Z must be >= 0" ); 129 endif 130 131 ## Initialize results 132 p = zeros( K, N+1 ); # p(k,i+1) is the probability that there are i jobs at server k, given that the network population is j 133 p(:,1) = 1; 134 U = R = Q = X = zeros( 1, K ); 135 X_s = 0; # System throughput 136 137 ## Main MVA loop, iterates over the population size 138 for n=1:N # over population size 139 140 j=1:n; 141 ## for i=1:K 142 ## R(i) = sum( j.*S(i,j).*p(i,j) ); 143 ## endfor 144 R = sum( repmat(j,K,1).*S(:,1:n).*p(:,1:n), 2)'; 145 146 R_s = dot( V, R ); # System response time 147 X_s = n / (Z+R_s); # System Throughput 148 ## G_N = G_Nm1 / X_s; G_Nm1 = G_N; 149 150 ## prepare for next iteration 151 for i=1:K 152 p(i, 1+j) = X_s * S(i,j) .* p(i,j) * V(i); 153 p(i, 1) = 1-sum(p(i,1+j)); 154 endfor 155 endfor 156 Q = X_s * ( V .* R ); 157 U = 1-p(:,1)'; # Service centers utilization 158 X = X_s * V; # Service centers throughput 159endfunction 160 161## Check degenerate case of N==0 (general LD case) 162%!test 163%! N = 0; 164%! S = [1 2; 3 4; 5 6; 7 8]; 165%! V = [1 1 1 4]; 166%! [U R Q X] = qncsmvald(N, S, V); 167%! assert( U, 0*V ); 168%! assert( R, 0*V ); 169%! assert( Q, 0*V ); 170%! assert( X, 0*V ); 171 172%!test 173%! # Exsample 3.42 p. 577 Jain 174%! V = [ 16 10 5 ]; 175%! N = 20; 176%! S = [ 0.125 0.3 0.2 ]; 177%! Sld = repmat( S', 1, N ); 178%! Z = 4; 179%! [U1 R1 Q1 X1] = qncsmvald(N,Sld,V,Z); 180%! [U2 R2 Q2 X2] = qncsmva(N,S,V,ones(1,3),Z); 181%! assert( U1, U2, 1e-3 ); 182%! assert( R1, R2, 1e-3 ); 183%! assert( Q1, Q2, 1e-3 ); 184%! assert( X1, X2, 1e-3 ); 185 186%!test 187%! # Example 8.7 p. 349 Bolch et al. 188%! N = 3; 189%! Sld = 1 ./ [ 2 4 4; ... 190%! 1.667 1.667 1.667; ... 191%! 1.25 1.25 1.25; ... 192%! 1 2 3 ]; 193%! V = [ 1 .5 .5 1 ]; 194%! [U R Q X] = qncsmvald(N,Sld,V); 195%! assert( Q, [0.624 0.473 0.686 1.217], 1e-3 ); 196%! assert( R, [0.512 0.776 1.127 1], 1e-3 ); 197%! assert( all( U<=1) ); 198 199%!test 200%! # Example 8.4 p. 333 Bolch et al. 201%! N = 3; 202%! S = [ .5 .6 .8 1 ]; 203%! m = [2 1 1 N]; 204%! Sld = zeros(3,N); 205%! Sld(1,:) = .5 ./ [1 2 2]; 206%! Sld(2,:) = [.6 .6 .6]; 207%! Sld(3,:) = [.8 .8 .8]; 208%! Sld(4,:) = 1 ./ [1 2 3]; 209%! V = [ 1 .5 .5 1 ]; 210%! [U1 R1 Q1 X1] = qncsmvald(N,Sld,V); 211%! [U2 R2 Q2 X2] = qncsmva(N,S,V,m); 212%! ## Note that qncsmvald computes the utilization in a different 213%! ## way as qncsmva; in fact, qncsmva knows that service 214%! ## center i has m(i)>1 servers, but qncsmvald does not. Thus, 215%! ## utilizations for multiple-server nodes cannot be compared 216%! assert( U1([2,3]), U2([2,3]), 1e-3 ); 217%! assert( R1, R2, 1e-3 ); 218%! assert( Q1, Q2, 1e-3 ); 219%! assert( X1, X2, 1e-3 ); 220 221