SMCTP2RFM Calculates the rainflow matrix for a SMCTP. CALL: [Frfc,mu_rfc] = smctp2rfm(P,F,c_m); Frfc = Rainflow Matrix [n,n] mu_rfc = Rainflow Counting intensity [n,n] P = Transition matrix for regime process. [r,r] F = Cell array of min-Max and Max-min matrices {r,2} F{i,1} = min-Max matrix, process i [n,n] F{i,2} = Max-min matrix, process i [n,n] c_m = Intensity of local minima, switching proc. [1,1] (Default: 1) If a matrix F{i,2}=[], then the process will be assumed to be time-reversible. Calculates the rainflow matrix for a switching process with a Markov chain of turning points within each regime. Example: (Example 4.1 in PhD thesis) P = [0.9 0.1; 0.05 0.95]; param = [-1 1 32]; u = levels(param); F1 = mktestmat(param,[-0.4 -0.3],0.15,1); F2 = mktestmat(param,[0.3 0.4],0.15,1); Frfc = smctp2rfm(P,{F1 F1'; F2 F2'}); cmatplot(u,u,Frfc) See also mctp2rfm, smctp2arfm
Calculates a counting distribution from a cycle matrix. | |
Calculates the stationary distribution for a Markov chain. | |
Current date and time as date vector. | |
Display message and abort function. | |
Matrix inverse. | |
Kronecker tensor product. | |
Extract lower triangular part. | |
Extract upper triangular part. | |
Display warning message; disable or enable warning messages. |
Auxiliary function used by ESTSMCTP | |
Script to computer exercises 2 | |
Script to computer exercises 4 | |
Calculates the rainflow matrix for a MCTP. | |
Rainflow matrix for Switching Markov Chains of Turning Points. | |
Quick test of the routines in module 'markov' |
001 function [F_rfc,mu_rfc,T] = smctp2rfm(P,F,c_m) 002 % SMCTP2RFM Calculates the rainflow matrix for a SMCTP. 003 % 004 % CALL: [Frfc,mu_rfc] = smctp2rfm(P,F,c_m); 005 % 006 % Frfc = Rainflow Matrix [n,n] 007 % mu_rfc = Rainflow Counting intensity [n,n] 008 % 009 % P = Transition matrix for regime process. [r,r] 010 % F = Cell array of min-Max and Max-min matrices {r,2} 011 % F{i,1} = min-Max matrix, process i [n,n] 012 % F{i,2} = Max-min matrix, process i [n,n] 013 % c_m = Intensity of local minima, switching proc. [1,1] 014 % (Default: 1) 015 % If a matrix F{i,2}=[], then the process will 016 % be assumed to be time-reversible. 017 % 018 % Calculates the rainflow matrix for a switching process 019 % with a Markov chain of turning points within each regime. 020 % 021 % Example: (Example 4.1 in PhD thesis) 022 % P = [0.9 0.1; 0.05 0.95]; 023 % param = [-1 1 32]; u = levels(param); 024 % F1 = mktestmat(param,[-0.4 -0.3],0.15,1); 025 % F2 = mktestmat(param,[0.3 0.4],0.15,1); 026 % Frfc = smctp2rfm(P,{F1 F1'; F2 F2'}); 027 % cmatplot(u,u,Frfc) 028 % 029 % See also mctp2rfm, smctp2arfm 030 031 % References 032 % 033 % P. Johannesson (1999): 034 % Rainflow Analysis of Switching Markov Loads. 035 % PhD thesis, Mathematical Statistics, Centre for Mathematical Sciences, 036 % Lund Institute of Technology. 037 % 038 % P. Johannesson (1998): 039 % Rainflow Cycles for Switching Processes with Markov Structure. 040 % Probability in the Engineering and Informational Sciences, 041 % Vol. 12, No. 2, pp. 143-175. 042 043 % Tested on Matlab 5.3 044 % 045 % History: 046 % Revised by PJ 23-Nov-1999 047 % updated for WAFO 048 % Created by PJ (Pär Johannesson) 1997 049 % Copyright (c) 1997-1998 by Pär Johannesson 050 % Toolbox: Rainflow Cycles for Switching Processes V.1.1, 22-Jan-1998 051 052 053 % Check input arguments 054 055 ni = nargin; 056 no = nargout; 057 error(nargchk(2,3,ni)); 058 059 if ni < 3 060 c_m=[]; 061 end 062 063 if isempty(c_m) 064 c_m=1; 065 end 066 067 % Define 068 069 Zstr = '123456789'; Zstr=Zstr'; 070 071 T(1,:)=clock; 072 073 r = length(P); % Number of regime states 074 075 n = length(F{1,1}); % Number of levels 076 077 % Check that the rowsums of P are equal to 1 078 079 sumP = sum(P'); 080 if sum(sumP == 1) ~= r 081 warning(['Rowsums of P not equal to 1. Renormalizing!']); 082 for i = 1:length(P) 083 P(i,:) = P(i,:)/sumP(i); 084 end 085 end 086 087 T(2,:)=clock; 088 089 % Normalize the rowsums of F{1,1},...,F{r,1} to 1 090 % ==> Q{1,1},...,Q{r,1} 091 092 for i = 1:r 093 QQ{i,1} = triu(F{i,1},1); % Set zeros below diagonal and diagonal 094 % Set negative elements to zero 095 [I,J] = find(QQ{i,1}<0); 096 if length(I) ~= 0 097 warning(['Negative elements in Q' Zstr(i) '. Setting to zero!']); 098 for k = 1:length(I) 099 QQ{i,1}(I(k),J(k)) = 0; 100 end 101 end 102 103 sumQQi = sum(QQ{i,1}'); 104 % Normalize rowsums 105 if sum(sumQQi == 1) ~= length(QQ{i,1}) 106 %disp(['Warning: Rowsums of Q' Zstr(i) ' not equal to 1. Renormalizing!']); 107 for j = 1:n-1 108 if sumQQi(j)~=0, QQ{i,1}(j,:) = QQ{i,1}(j,:)/sumQQi(j); end 109 end 110 end 111 end 112 113 T(3,:)=clock; 114 115 % Normalize the rowsums of F{1,2},...,F{r,2} to 1 116 % ==> Q{1,2},...,Q{r,2} 117 118 % Normalize the rowsums of Fh1,...,Fhr to 1 ==> Q1,...,Qr 119 120 for i = 1:r 121 if isempty(F{i,2}) % Time-reversible 122 QQ{i,2} = F{i,1}'; 123 else % Fhi is given 124 QQ{i,2} = F{i,2}; 125 end 126 127 QQ{i,2} = tril(QQ{i,2},-1); % Set zeros above diagonal and diagonal 128 % Set negative elements to zero 129 [I,J] = find(QQ{i,2}<0); 130 if length(I) ~= 0 131 warning(['Negative elements in Qh' Zstr(i) '. Setting to zero!']); 132 for k = 1:length(I) 133 QQ{i,2}(I(k),J(k)) = 0; 134 end 135 end 136 137 sumQQi = sum(QQ{i,2}'); 138 if sum(sumQQi == 1) ~= length(QQ{i,2}) 139 %disp(['Warning: Rowsums of Qh' Zstr(i) ' not equal to 1. Renormalizing!']); 140 for j = 2:n 141 if sumQQi(j)~=0, QQ{i,2}(j,:) = QQ{i,2}(j,:)/sumQQi(j); end 142 end 143 end 144 end 145 146 T(4,:)=clock; 147 148 % Make the transition matrix Q for the joint min-Max process 149 150 Q = zeros(n*r,n*r); 151 I = 0:r:(n-1)*r; 152 for z=1:r 153 Q0 = kron(QQ{z,1},P); 154 Q(I+z,:) = Q0(I+z,:); 155 end 156 157 T(5,:)=clock; 158 159 % Make the transition matrix Qh for the joint Max-min process 160 161 Qh = zeros(n*r,n*r); 162 I = 0:r:(n-1)*r; 163 for z=1:r 164 Q0 = kron(QQ{z,2},P); 165 Qh(I+z,:) = Q0(I+z,:); 166 end 167 168 T(6,:)=clock; 169 170 % Stationary distribution (=ro) of local minima with transition matrix 171 % Qt = Q*Qh = "Transition matrix for min-to-min" 172 173 Qt = Q*Qh; 174 ro = mc2stat(Qt(1:r*(n-1),1:r*(n-1))); % Stationary distr., row vector 175 ro = [ro zeros(1,r)]; % Minimum can't reach the highest level 176 177 T(7,:)=clock; 178 179 % Calculate rainflow matrix 180 181 mu_rfc = zeros(n,n); 182 EYE = eye(n*r,n*r); 183 Qcumsum = fliplr(cumsum(fliplr(Q)')'); 184 185 %fprintf(1,'Calculating row '); 186 for i=2:n 187 %fprintf(1,'-%1d',i); 188 for j=i-1:min([i n-1]) 189 q = Qcumsum(:,r*j+1); 190 191 d = q(1:r*(i-1)); % 1:i-1 192 Ro = ro(1:r*(i-1)); % 1:i-1 193 194 mu_rfc(i,j) = Ro*d; 195 end 196 for j=i+1:n-1 197 q = Qcumsum(:,r*j+1); 198 199 A = Q(r*(i-1)+1:r*(j-1),r*(i+1-1)+1:r*j); % i:j-1, i+1:j 200 Ah = Qh(r*(i+1-1)+1:r*j,r*(i-1)+1:r*(j-1)); % i+1:j, i:j-1 201 C = Q(1:r*(i-1),r*(i+1-1)+1:r*j); % 1:i-1, i+1:j 202 Eye = EYE(1:r*(j-i),1:r*(j-i)); % eye (size of A*Ah) 203 d = q(1:r*(i-1)); % 1:i-1 204 e = q(r*(i-1)+1:r*(j-1)); % i:j-1 205 Ro = ro(1:r*(i-1)); % 1:i-1 206 207 mu_rfc(i,j) = Ro*(d+C*Ah*inv(Eye-A*Ah)*e); 208 end 209 end 210 %fprintf(1,'\n'); 211 212 213 % Calculation of the intensity of local minima for the swithching process 214 215 mu_rfc = c_m*mu_rfc; 216 217 T(8,:)=clock; 218 219 % Convert: mu_rfc --> F_rfc 220 221 F_rfc = zeros(n,n); 222 for i= 1:n-1 223 for j= i+1:n 224 F_rfc(i,j) = mu_rfc(i+1,j-1) - mu_rfc(i,j-1) - mu_rfc(i+1,j) + mu_rfc(i,j); 225 end 226 end 227 228 I = F_rfc<0; 229 if sum(sum(I)) ~= 0 230 warning(['Negative elements in calculated rainflow matrix F_rfc. Setting to zero!']); 231 F_rfc(I) = 0; 232 end 233 234 235 if no>1 236 mu_rfc = cmat2nt(F_rfc); % Calculate rainflow counting intensity 237 end 238 239 240 T(9,:)=clock; 241 242
Comments or corrections to the WAFO group