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