SMC2RFM Calculates the rainflow matrix/intensity for a switching Markov chain. [F_rfc,mu_rfc] = smc2rfc(P,Q,1); [F_rfc,mu_rfc] = smc2rfc(P,Q,[2 h]); F_rfc = Rainflow matrix / Rainflow intensity [NxN] mu_rfc = Rainflow counting intensity [NxN] P = Transition matrix for regime process [rxr] Q = Cell array of transition matrices {r,1} Q{i} = Transition matrix for Markov chain i [nxn] def = Definition 1: Markov chain (default), N=n 2: Discretized Markov chain, N=n+1 h = Discretization step (ONLY Def 2!) Calculates (1) the rainflow matrix for a switching Markov chain OR (2) the rainflow intensity for a discretized switching Markov chain. Example: param = [-1 1 32]; u = levels(param); [F1,F2] = mktestmat(param,[-0.3 0.3],0.15,1,-Inf); Frfc = smc2rfm(P,{F1;F2}); cmatplot(u,u,Frfc) See also mc2rfm, smctp2rfm, smc2stat, cmatplot
Calculates a counting distribution from a cycle matrix. | |
Calculates the stationary distribution for a Markov chain. | |
Display message and abort function. | |
Matrix inverse. | |
Kronecker tensor product. | |
Display warning message; disable or enable warning messages. |
Calculates the rainflow matrix/intensity for a Markov chain. | |
Demo for switching AR(1)-processes. | |
Quick test of the routines in module 'markov' |
001 function [F_rfc,mu_rfc] = smc2rfc(P,Qc,def) 002 %SMC2RFM Calculates the rainflow matrix/intensity for a switching Markov chain. 003 % 004 % [F_rfc,mu_rfc] = smc2rfc(P,Q,1); 005 % [F_rfc,mu_rfc] = smc2rfc(P,Q,[2 h]); 006 % 007 % F_rfc = Rainflow matrix / Rainflow intensity [NxN] 008 % mu_rfc = Rainflow counting intensity [NxN] 009 % 010 % P = Transition matrix for regime process [rxr] 011 % Q = Cell array of transition matrices {r,1} 012 % Q{i} = Transition matrix for Markov chain i [nxn] 013 % def = Definition 1: Markov chain (default), N=n 014 % 2: Discretized Markov chain, N=n+1 015 % h = Discretization step (ONLY Def 2!) 016 % 017 % Calculates 018 % (1) the rainflow matrix for a switching Markov chain OR 019 % (2) the rainflow intensity for a discretized switching Markov chain. 020 % 021 % Example: 022 % param = [-1 1 32]; u = levels(param); 023 % [F1,F2] = mktestmat(param,[-0.3 0.3],0.15,1,-Inf); 024 % Frfc = smc2rfm(P,{F1;F2}); 025 % cmatplot(u,u,Frfc) 026 % 027 % See also mc2rfm, smctp2rfm, smc2stat, cmatplot 028 029 % References 030 % 031 % P. Johannesson (1999): 032 % Rainflow Analysis of Switching Markov Loads. 033 % PhD thesis, Mathematical Statistics, Centre for Mathematical Sciences, 034 % Lund Institute of Technology. 035 % 036 % P. Johannesson (1998): 037 % Rainflow Cycles for Switching Processes with Markov Structure. 038 % Probability in the Engineering and Informational Sciences, 039 % Vol. 12, No. 2, pp. 143-175. 040 041 % Tested on Matlab 5.3 042 % 043 % History: 044 % Revised by PJ 19-May-2000 045 % Changer disp(...) to warning(...) . 046 % Revised by PJ 23-Nov-1999 047 % updated for WAFO 048 % Created by PJ (Pär Johannesson) 1997 049 % Copyright (c) 1997 by Pär Johannesson 050 % Toolbox: Rainflow Cycles for Switching Processes V.1.0, 2-Oct-1997 051 052 053 % Check input arguments 054 055 ni = nargin; 056 no = nargout; 057 error(nargchk(2,3,ni)); 058 059 if ni<3, def = []; end 060 if isempty(def), def = 1; end 061 062 % Define 063 064 Zstr = '123456789'; 065 066 r = length(P); % Number of regime states 067 n = length(Qc{1}); % Number of levels 068 069 % Check that the rowsums of P are equal to 1 070 071 sumP = sum(P'); 072 if sum(sumP == 1) ~= length(P) 073 warning(['Rowsums of P not equal to 1. Renormalizing.']); 074 for i = 1:length(P) 075 P(i,:) = P(i,:)/sumP(i); 076 end 077 end 078 079 % Check that the rowsums of Qc{1},...,Qc{r} are equal to 1 080 081 for i = 1:r 082 sumQi = sum(Qc{i}'); 083 if sum(sumQi == 1) ~= length(Qc{i}) 084 warning(['Rowsums of Q{' Zstr(i) '} not equal to 1. Renormalizing.']); 085 for j = 1:length(Qc{i}) 086 Qc{i}(j,:) = Qc{i}(j,:)/sumQi(j); 087 end 088 end 089 end 090 091 092 % Make the transition matrix Q for the joint MC (X_k,Z_k) 093 094 Q = zeros(n*r,n*r); 095 I = 0:r:(n-1)*r; 096 for z=1:r 097 QQ = kron(Qc{z},P); 098 Q(I+z,:) = QQ(I+z,:); 099 end 100 101 102 % Stationary distribution = ro of Q 103 104 ro = mc2stat(Q); 105 106 % Calculate rainflow matrix / rainflow intensity 107 108 if def(1) == 1 % Markov Chain 109 110 N = n; 111 mu_rfc = zeros(N,N); 112 EYE = eye(n*r,n*r); 113 Qcumsum = fliplr(cumsum(fliplr(Q)')'); 114 115 for i=2:n 116 for j=i-1:n-1 117 q = Qcumsum(:,r*j+1); 118 119 A = Q(r*(i-1)+1:r*j,r*(i-1)+1:r*j); % i:j, i:j 120 C = Q(1:r*(i-1),r*(i-1)+1:r*j); % 1:i-1, i:j 121 Eye = EYE(1:r*(j-i+1),1:r*(j-i+1)); % eye (size of A) 122 d = q(1:r*(i-1)); % 1:i-1 123 e = q(r*(i-1)+1:r*j); % i:j 124 Ro = ro(1:r*(i-1)); % 1:i-1 125 126 if j == i-1 % isempty(A) 127 mu_rfc(i,j) = Ro*d; 128 else 129 mu_rfc(i,j) = Ro*(d+C*inv(Eye-A)*e); 130 end 131 end 132 end 133 134 elseif def(1) == 2 % Discretized Markov Process 135 136 N = n+1; 137 mu_rfc = zeros(N,N); 138 EYE = eye(n*r,n*r); 139 Qcumsum = fliplr(cumsum(fliplr(Q)')'); 140 141 for i=2:N-1 142 for j=i:N-1 143 q = Qcumsum(:,r*(j-1)+1); 144 145 A = Q(r*(i-1)+1:r*(j-1),r*(i-1)+1:r*(j-1)); % i:j-1, i:j-1 146 C = Q(1:r*(i-1),r*(i-1)+1:r*(j-1)); % 1:i-1, i:j-1 147 Eye = EYE(1:r*(j-i),1:r*(j-i)); % eye (size of A) 148 d = q(1:r*(i-1)); % 1:i-1 149 e = q(r*(i-1)+1:r*(j-1)); % i-1:j-1 150 Ro = ro(1:r*(i-1)); % 1:i-1 151 152 if i == j % isempty(A) 153 mu_rfc(i,j) = Ro*d; 154 else 155 mu_rfc(i,j) = Ro*(d+C*inv(Eye-A)*e); 156 end 157 end 158 end 159 160 end 161 162 % Fill the missing triangel in the mu_rfc matrix 163 if def(1) == 1 164 165 F_rfc = zeros(n); 166 for i = 1:n-1 167 for j= i+1:n 168 F_rfc(i,j) = mu_rfc(i+1,j-1)-mu_rfc(i,j-1)-mu_rfc(i+1,j)+mu_rfc(i,j); 169 end 170 end 171 172 mu_rfc = cmat2nt(F_rfc); 173 174 % Old version 175 % mu_rfc=flipud(mu_rfc'); % Convert to WAT matrix-format 176 % for k = 2:N-1 177 % for m = 1:N-k 178 % i = k+m; 179 % j = N+1+k-i; 180 % mu_rfc(i,j) = mu_rfc(i,j-1)+mu_rfc(i-1,j)-mu_rfc(i-1,j-1); 181 % end 182 % end 183 % 184 % mu_rfc = flipud(mu_rfc)'; % Convert to PJ-def 185 % 186 % F_rfc = nt2cmat(mu_rfc); % Calculate Rainflow matrix 187 % 188 %% F_rfc = fliplr(triu(fliplr(F_rfc),1)); % Zeros blow SW-NE diagonal 189 190 end 191 192 if def(1) == 2 193 194 mu_rfc=flipud(mu_rfc'); % Convert to WAT matrix-format 195 for k = 1:N-1 196 for m = 1:N-k 197 i = k+m; 198 j = N+1+k-i; 199 mu_rfc(i,j) = mu_rfc(i,j-1)+mu_rfc(i-1,j)-mu_rfc(i-1,j-1); 200 end 201 end 202 203 % Calculate rainflow intensity by numerical differentiation 204 205 h = def(2); 206 F_rfc = zeros(N,N); 207 208 for i=2:N-1 209 for j=2:N-i 210 F_rfc(i,j) = (mu_rfc(i-1,j-1)-mu_rfc(i-1,j+1) -(mu_rfc(i+1,j-1)-mu_rfc(i+1,j+1))) / (2*h)^2; 211 end 212 end 213 214 F_rfc = flipud(F_rfc)'; % Convert to PJ-def 215 mu_rfc = flipud(mu_rfc)'; % Convert to PJ-def 216 217 end 218 219 220 221
Comments or corrections to the WAFO group