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{L} =} ctmcexps (@var{Q}, @var{t}, @var{p} ) 21## @deftypefnx {Function File} {@var{L} =} ctmcexps (@var{Q}, @var{p}) 22## 23## @cindex Markov chain, continuous time 24## @cindex continuous time Markov chain 25## @cindex expected sojourn time, CTMC 26## @cindex CTMC 27## 28## With three arguments, compute the expected times @code{@var{L}(i)} 29## spent in each state @math{i} during the time interval @math{[0,t]}, 30## assuming that the initial occupancy vector is @var{p}. With two 31## arguments, compute the expected time @code{@var{L}(i)} spent in each 32## transient state @math{i} until absorption. 33## 34## @strong{Note:} In its current implementation, this function 35## requires that an absorbing state is reachable from any 36## non-absorbing state of @math{Q}. 37## 38## @strong{INPUTS} 39## 40## @table @code 41## 42## @item @var{Q}(i,j) 43## @math{N \times N} infinitesimal generator matrix. @code{@var{Q}(i,j)} 44## is the transition rate from state @math{i} to state @math{j}, 45## @math{1 @leq{} i, j @leq{} N}, @math{i \neq j}. 46## The matrix @var{Q} must also satisfy the 47## condition @math{\sum_{j=1}^N Q_{i,j} = 0} for every @math{i=1, @dots{}, N}. 48## 49## @item @var{t} 50## If given, compute the expected sojourn times in @math{[0,t]} 51## 52## @item @var{p}(i) 53## Initial occupancy probability vector; @code{@var{p}(i)} is the 54## probability the system is in state @math{i} at time 0, @math{i = 1, 55## @dots{}, N} 56## 57## @end table 58## 59## @strong{OUTPUTS} 60## 61## @table @code 62## 63## @item @var{L}(i) 64## If this function is called with three arguments, @code{@var{L}(i)} is 65## the expected time spent in state @math{i} during the interval 66## @math{[0,t]}. If this function is called with two arguments 67## @code{@var{L}(i)} is the expected time spent in transient state 68## @math{i} until absorption; if state @math{i} is absorbing, 69## @code{@var{L}(i)} is zero. 70## 71## @end table 72## 73## @seealso{dtmcexps} 74## 75## @end deftypefn 76 77## Author: Moreno Marzolla <moreno.marzolla(at)unibo.it> 78## Web: http://www.moreno.marzolla.name/ 79 80function L = ctmcexps( Q, varargin ) 81 82 persistent epsilon = 10*eps; 83 84 if ( nargin < 2 || nargin > 3 ) 85 print_usage(); 86 endif 87 88 [N err] = ctmcchkQ(Q); 89 90 (N>0) || ... 91 error(err); 92 93 if ( nargin == 2 ) 94 p = varargin{1}; 95 else 96 t = varargin{1}; 97 p = varargin{2}; 98 endif 99 100 ( isvector(p) && length(p) == size(Q,1) && all(p>=0) && abs(sum(p)-1.0)<epsilon ) || ... 101 error( "p must be a probability vector" ); 102 103 p = p(:)'; # make p a row vector 104 105 if ( nargin == 3 ) # non-absorbing case 106 if ( isscalar(t) ) 107 (t >= 0 ) || ... 108 error( "t must be >= 0" ); 109 ## F(x) are the transient state occupancy probabilities at time x 110 ## F(x) = p*expm(Q*x) (see function ctmc()). 111 F = @(x) (p*expm(Q*x)); 112 L = quadv(F,0,t); 113 else 114 ## FIXME: deprecate this? 115 ( isvector(t) && abs(t(1)) < epsilon ) || ... 116 error( "t must be a vector, and t(1) must be 0.0" ); 117 t = t(:)'; # make t a row vector 118 ff = @(x,t) (x(:)'*Q+p); 119 fj = @(x,t) (Q); 120 L = lsode( {ff, fj}, zeros(size(p)), t ); 121 endif 122 else # absorbing case 123 124 ## Identify absorbing states. If there are no absorbing states, 125 ## raise an error. 126 127 N = rows(Q); 128 tr = find( any( abs(Q) > epsilon, 2 ) ); # non-absorbing states 129 if ( length( tr ) == N ) 130 error( "There are no absorbing states" ); 131 endif 132 133 QN = Q(tr,tr); 134 pN = p(tr); 135 LN = -pN*inv(QN); 136 L = zeros(1,N); 137 L(tr) = LN; 138 endif 139endfunction 140%!test 141%! Q = [-1 1; 1 -1]; 142%! L = ctmcexps(Q,10,[1 0]); 143%! L = ctmcexps(Q,linspace(0,10,100),[1 0]); 144 145%!test 146%! Q = ctmcbd( [1 2 3], [3 2 1] ); 147%! p0 = [1 0 0 0]; 148%! t = linspace(0,10,10); 149%! L1 = L2 = zeros(length(t),4); 150%! # compute L using the differential equation formulation 151%! ff = @(x,t) (x(:)'*Q+p0); 152%! fj = @(x,t) (Q); 153%! L1 = lsode( {ff, fj}, zeros(size(p0)), t ); 154%! # compute L using ctmcexps (integral formulation) 155%! for i=1:length(t) 156%! L2(i,:) = ctmcexps(Q,t(i),p0); 157%! endfor 158%! assert( L1, L2, 1e-5); 159 160%!demo 161%! lambda = 0.5; 162%! N = 4; 163%! b = lambda*[1:N-1]; 164%! d = zeros(size(b)); 165%! Q = ctmcbd(b,d); 166%! t = linspace(0,10,100); 167%! p0 = zeros(1,N); p0(1)=1; 168%! L = zeros(length(t),N); 169%! for i=1:length(t) 170%! L(i,:) = ctmcexps(Q,t(i),p0); 171%! endfor 172%! plot( t, L(:,1), ";State 1;", "linewidth", 2, ... 173%! t, L(:,2), ";State 2;", "linewidth", 2, ... 174%! t, L(:,3), ";State 3;", "linewidth", 2, ... 175%! t, L(:,4), ";State 4;", "linewidth", 2 ); 176%! legend("location","northwest"); legend("boxoff"); 177%! xlabel("Time"); 178%! ylabel("Expected sojourn time"); 179 180%!demo 181%! lambda = 0.5; 182%! N = 4; 183%! b = lambda*[1:N-1]; 184%! d = zeros(size(b)); 185%! Q = ctmcbd(b,d); 186%! p0 = zeros(1,N); p0(1)=1; 187%! L = ctmcexps(Q,p0); 188%! disp(L); 189