1## Copyright (C) 2008, 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{t} =} ctmcmtta (@var{Q}, @var{p}) 21## 22## @cindex Markov chain, continuous time 23## @cindex continuous time Markov chain 24## @cindex CTMC 25## @cindex mean time to absorption, CTMC 26## 27## Compute the Mean-Time to Absorption (MTTA) of the CTMC described by 28## the infinitesimal generator matrix @var{Q}, starting from initial 29## occupancy probabilities @var{p}. If there are no absorbing states, this 30## function fails with an error. 31## 32## @strong{INPUTS} 33## 34## @table @code 35## 36## @item @var{Q}(i,j) 37## @math{N \times N} infinitesimal generator matrix. @code{@var{Q}(i,j)} 38## is the transition rate from state @math{i} to state @math{j}, @math{i 39## \neq j}. The matrix @var{Q} must satisfy the condition 40## @math{\sum_{j=1}^N Q_{i,j} = 0} 41## 42## @item @var{p}(i) 43## probability that the system is in state @math{i} 44## at time 0, for each @math{i=1, @dots{}, N} 45## 46## @end table 47## 48## @strong{OUTPUTS} 49## 50## @table @code 51## 52## @item @var{t} 53## Mean time to absorption of the process represented by matrix @var{Q}. 54## If there are no absorbing states, this function fails. 55## 56## @end table 57## 58## @strong{REFERENCES} 59## 60## @itemize 61## @item 62## G. Bolch, S. Greiner, H. de Meer and 63## K. Trivedi, @cite{Queueing Networks and Markov Chains: Modeling and 64## Performance Evaluation with Computer Science Applications}, Wiley, 65## 1998. 66## @end itemize 67## 68## @seealso{ctmcexps} 69## 70## @end deftypefn 71 72## Author: Moreno Marzolla <moreno.marzolla(at)unibo.it> 73## Web: http://www.moreno.marzolla.name/ 74 75function t = ctmcmtta( Q, p ) 76 77 persistent epsilon = 10*eps; 78 79 if ( nargin != 2 ) 80 print_usage(); 81 endif 82 83 [N err] = ctmcchkQ(Q); 84 85 (N>0) || ... 86 error(err); 87 88 ( isvector(p) && length(p) == N && all(p>=0) && abs(sum(p)-1.0)<epsilon ) || ... 89 error( "p must be a probability vector" ); 90 p = p(:)'; 91 92 L = ctmcexps(Q,p); 93 t = sum(L); 94endfunction 95%!test 96%! Q = [0 1 0; 1 0 1; 0 1 0 ]; Q -= diag( sum(Q,2) ); 97%! fail( "ctmcmtta(Q,[1 0 0])", "no absorbing"); 98 99%!test 100%! Q = [0 1 0; 1 0 1; 0 0 0; 0 0 0 ]; 101%! fail( "ctmcmtta(Q,[1 0 0])", "square matrix"); 102 103%!test 104%! Q = [0 1 0; 1 0 1; 0 0 0 ]; 105%! fail( "ctmcmtta(Q,[1 0 0])", "infinitesimal"); 106 107%!test 108%! Q = [ 0 0.1 0 0; ... 109%! 0.9 0 0.1 0; ... 110%! 0 0.9 0 0.1; ... 111%! 0 0 0 0 ]; 112%! Q -= diag( sum(Q,2) ); 113%! assert( ctmcmtta( Q,[0 0 0 1] ), 0 ); # state 4 is absorbing 114 115%!test 116%! Q = [-1 1; 0 0]; 117%! assert( ctmcmtta( Q, [0 1] ), 0 ); # state 2 is absorbing 118%! assert( ctmcmtta( Q, [1 0] ), 1 ); # the result has been computed by hand 119 120## Compute the MTTA of a pure death process with 4 states 121## (state 1 is absorbing). State 4 is the initial state. 122%!demo 123%! mu = 0.01; 124%! death = [ 3 4 5 ] * mu; 125%! birth = 0*death; 126%! Q = ctmcbd(birth,death); 127%! t = ctmcmtta(Q,[0 0 0 1]) 128 129%!demo 130%! N = 100; 131%! birth = death = ones(1,N-1); birth(1) = death(N-1) = 0; 132%! Q = diag(birth,1)+diag(death,-1); 133%! Q -= diag(sum(Q,2)); 134%! t = zeros(1,N/2); 135%! initial_state = 1:(N/2); 136%! for i=initial_state 137%! p = zeros(1,N); p(i) = 1; 138%! t(i) = ctmcmtta(Q,p); 139%! endfor 140%! plot(initial_state,t,"+"); 141%! xlabel("Initial state"); 142%! ylabel("MTTA"); 143 144 145