MCTP2ARFM Calculates asymmetric rainflow matrix for a MCTP. CALL: [F_rfc] = mctp2arfm(F); [F_rfc,FF,FFr] = mctp2arfm(F,c_m); F_rfc = Asymmetric Rainflow Matrix [n,n] FF = Cell array of min-Max and Max-min matrices {1,2} FFr = Cell array of min-Max and Max-min matrices {1,2} F = Cell array of min-Max and Max-min matrices {1,2} F{1,1} = min-Max matrix [n,n] F{1,2} = Max-min matrix [n,n] c_m = Intensity of local minima, switching proc. [1,1] (Default: 1) Calculates the expected Asymmetric RainFlow Matrix (ARFM) for a Markov Chain of Turning Points. Recursive formulation a'la Igor If a matrix F{1,2}=[], then the process will be assumed to be time-reversible. Example: param = [-1 1 32]; u = levels(param); F = mktestmat(param,[-0.2 0.2],0.15,2); Frfc = mctp2rfm({F []}); Farfc = mctp2arfm({F []}); cmatplot(u,u,{Frfc Farfc},3) sum(sum(abs((Farfc+Farfc')-(Frfc+Frfc')))) % should be zero See also arfm2mctp, smctp2arfm, mctp2rfm, cmatplot
Converts a matrix to a transition matrix. | |
Calculates the stationary distribution for a Markov chain. | |
Current date and time as date vector. | |
Display message and abort function. | |
Matrix inverse. | |
Display warning message; disable or enable warning messages. |
Quick test of the routines in module 'markov' |
0001 function [F_rfc,FF,FFr,T] = mctp2arfm(F,c_m) 0002 %MCTP2ARFM Calculates asymmetric rainflow matrix for a MCTP. 0003 % 0004 % CALL: [F_rfc] = mctp2arfm(F); 0005 % [F_rfc,FF,FFr] = mctp2arfm(F,c_m); 0006 % 0007 % F_rfc = Asymmetric Rainflow Matrix [n,n] 0008 % FF = Cell array of min-Max and Max-min matrices {1,2} 0009 % FFr = Cell array of min-Max and Max-min matrices {1,2} 0010 % 0011 % F = Cell array of min-Max and Max-min matrices {1,2} 0012 % F{1,1} = min-Max matrix [n,n] 0013 % F{1,2} = Max-min matrix [n,n] 0014 % c_m = Intensity of local minima, switching proc. [1,1] 0015 % (Default: 1) 0016 % 0017 % Calculates the expected Asymmetric RainFlow Matrix (ARFM) for a 0018 % Markov Chain of Turning Points. Recursive formulation a'la Igor 0019 % If a matrix F{1,2}=[], then the process will be assumed to be 0020 % time-reversible. 0021 % 0022 % Example: 0023 % param = [-1 1 32]; u = levels(param); 0024 % F = mktestmat(param,[-0.2 0.2],0.15,2); 0025 % Frfc = mctp2rfm({F []}); 0026 % Farfc = mctp2arfm({F []}); 0027 % cmatplot(u,u,{Frfc Farfc},3) 0028 % sum(sum(abs((Farfc+Farfc')-(Frfc+Frfc')))) % should be zero 0029 % 0030 % See also arfm2mctp, smctp2arfm, mctp2rfm, cmatplot 0031 0032 % References: 0033 % 0034 % P. Johannesson (1999): 0035 % Rainflow Analysis of Switching Markov Loads. 0036 % PhD thesis, Mathematical Statistics, Centre for Mathematical Sciences, 0037 % Lund Institute of Technology. 0038 0039 % Tested on Matlab 5.3 0040 % 0041 % History: 0042 % Revised by PJ 09-Apr-2001 0043 % updated for WAFO 0044 % Created by PJ (Pär Johannesson) 1998 0045 % Copyright (c) 1997-1998 by Pär Johannesson 0046 % Toolbox: Rainflow Cycles for Switching Processes V.1.1, 22-Jan-1998 0047 0048 % This program used the formulation where the probabilities 0049 % of the events are calculated using "elementary" events for 0050 % the MCTP. 0051 % 0052 % Standing 0053 % pS = Max*pS1*pS2*pS3; 0054 % F_rfc(i,j) = pS; 0055 % Hanging 0056 % pH = Min*pH1*pH2*pH3; 0057 % F_rfc(j,i) = pH; 0058 % 0059 % The cond. prob. pS1, pS2, pS3, pH1, pH2, pH3 are calculated using 0060 % the elementary cond. prob. C, E, R, D, E3, Ch, Eh, Rh, Dh, E3h. 0061 0062 T(1,:)=clock; 0063 0064 % Check input arguments 0065 0066 ni = nargin; 0067 no = nargout; 0068 error(nargchk(1,2,ni)); 0069 0070 if ni < 2 0071 c_m=1; 0072 end 0073 0074 % Define 0075 0076 n = length(F{1,1}); % Number of levels 0077 0078 % Normalize the rowsums of F{1,1} to 1 ==> Q 0079 0080 Q = F{1,1}; 0081 Q = mat2tmat(Q,1); % Convert to min-max transition matrix 0082 0083 % Normalize the rowsums of F{1,2} to 1 ==> Qh 0084 0085 if isempty(F{1,2}) % Time-reversible 0086 Qh = F{1,1}'; 0087 else % Fh is given 0088 Qh = F{1,2}; 0089 end 0090 0091 Qh = mat2tmat(Qh,-1); % Convert to max-min transition matrix 0092 0093 T(2,:)=clock; 0094 0095 % Stationary distribution (=ro) of local minima with transition matrix 0096 % Qt = Q*Qh = "Transition matrix for min-to-min" 0097 0098 Qt = Q*Qh; 0099 ro = mc2stat(Qt(1:(n-1),1:(n-1))); % Stationary distr., row vector 0100 ro = [ro 0]; % Minimum can't reach the highest level 0101 I = find(ro<0); 0102 if ~isempty(I) 0103 warning(['Negative elements in ro. Setting to zero!']); 0104 ro(I) = zeros(size(I)); 0105 end 0106 0107 % Stationary distribution (=roh) of local maxima with transition matrix 0108 % Qt = Qh*Q = "Transition matrix for max-to-max" 0109 0110 Qth = Qh*Q; 0111 roh = mc2stat(Qth(2:n,2:n)); % Stationary distr., row vector 0112 roh = [0 roh]; % Maximum can't reach the highest level 0113 I = find(roh<0); 0114 if ~isempty(I) 0115 warning(['Negative elements in roh. Setting to zero!']); 0116 roh(I) = zeros(size(I)); 0117 end 0118 0119 T(3,:)=clock; 0120 0121 % Create Transition matrices for time-reversed MC 0122 0123 % Backward min-to-max 0124 I1 = find(ro>0); 0125 ro_inv = zeros(1,n); ro_inv(I1) = 1./ro(I1); 0126 Qr = Qh' .* (ro_inv'*roh); 0127 0128 % Backward max-to-min 0129 I1 = find(roh>0); 0130 roh_inv = zeros(1,n); roh_inv(I1) = 1./roh(I1); 0131 Qrh = Q' .* (roh_inv'*ro); 0132 0133 T(4,:)=clock; 0134 0135 % Define matrices for going to or above a level 0136 %R = fliplr(cumsum(fliplr(Q),2)); 0137 %Rr = fliplr(cumsum(fliplr(Qr),2)); 0138 0139 % Define matrices for going to or below a level 0140 %Rh = cumsum(Qh,2); 0141 %Rrh = cumsum(Qrh,2); 0142 0143 % Compure min-max and max-min matrices for model 0144 FF{1,1} = diag(ro)*Q; % min-max matrix 0145 FF{1,2} = diag(roh)*Qh; % max-min matrix 0146 FFr{1,1} = diag(roh)*Qrh; % min-max matrix (from reversed chain) 0147 FFr{1,2} = diag(ro)*Qr; % min-max matrix (from reversed chain) 0148 % FF and FFr should be identical 0149 0150 T(5,:)=clock; 0151 0152 % Calculate min-max and max-min matrices from input F 0153 %F1 = F{1,1}; 0154 %F1 = F1/sum(sum(F1)); 0155 %F1h = F{1,2}; 0156 %F1h = F1h/sum(sum(F1h)); 0157 % The model FF need not be equal to F 0158 F1=Q; F1h=Qh; 0159 0160 % Calculate rainflow matrix 0161 0162 F_rfc = zeros(n,n); 0163 EYE = eye(n,n); 0164 0165 %fprintf(1,'Calculating row '); 0166 for i=1:n-1 0167 % fprintf(1,'-%1d',i); 0168 0169 for j=i+1:n 0170 0171 xx = FF{1,1}(i,j); 0172 yy = FF{1,2}(j,i); 0173 0174 x = Q(i,j)*ro(i); 0175 y = Qh(j,i)*roh(j); 0176 0177 %fprintf(1,'x=%f, y=%f, xx=%f, yy=%f\n',x,y,xx,yy); 0178 0179 % Min = sum(F1(i,:)); 0180 % Max = sum(F1h(j,:)); 0181 Min = sum(F1(i,:)+F1h(:,i)')/2; 0182 Max = sum(F1h(j,:)+F1(:,j)')/2; 0183 0184 Ro = ro(i); % Probability of "min=i" 0185 Roh = roh(j); % Probability of "max=j" 0186 0187 % fprintf(1,'Min=%f, Max=%f, Ro=%f, Roh=%f\n',Min,Max,Ro,Roh); 0188 0189 Min = Ro; Max = Roh; % Just to be sure 0190 0191 if j-i == 1 % First subdiagonal 0192 0193 C = y/Min; Ch = x/Max; 0194 D = 0; Dh = 0; 0195 E = 1-y/Min; Eh = 1-x/Max; 0196 R = 0; Rh = 0; 0197 E3 = 1; E3h = 1; 0198 0199 elseif j-i == 2 % Second subdiagonal 0200 0201 % For Part 1 & 2 of cycle 0202 0203 I = i+1:j-2; 0204 J = i+2:j-1; 0205 0206 e = 1 - sum(Qr(I,i+2:j),2); 0207 eh = 1 - sum(Qrh(J,i:j-2),2); 0208 0209 C = y/Min; 0210 Ch = x/Max; 0211 d1 = Qr(i,i+1)*(1-Qrh(i+1,i)); 0212 D = d1; 0213 d1h = Qrh(j,j-1)*(1-Qr(j-1,j)); 0214 Dh = d1h; 0215 d0 = sum(Qr(i,i+1:j-1)); 0216 E = 1-d0-y/Min; 0217 d0h = sum(Qrh(j,i+1:j-1)); 0218 Eh = 1-d0h-x/Max; 0219 r1 = Qr(i,i+1)*Qrh(i+1,i); 0220 R = r1; 0221 r1h = Qrh(j,j-1)*Qr(j-1,j); 0222 Rh = r1h; 0223 0224 % For Part 3 of cycle 0225 0226 d3 = sum(Q(i,i+1:j-1)); 0227 E3 = 1-d3; 0228 d3h = sum(Qh(j,i+1:j-1)); 0229 E3h = 1-d3h; 0230 0231 else 0232 0233 Eye = EYE(1:j-i-2,1:j-i-2); 0234 0235 % For Part 1 & 2 of cycle 0236 0237 I = i+1:j-2; 0238 J = i+2:j-1; 0239 A = Qr(I,J); 0240 Ah = Qrh(J,I); 0241 a = Qr(i,J); 0242 ah = Qrh(j,I); 0243 b = Qr(I,j); 0244 bh = Qrh(J,i); 0245 0246 e = 1 - sum(Qr(I,i+2:j),2); 0247 eh = 1 - sum(Qrh(J,i:j-2),2); 0248 0249 Inv = inv(Eye-A*Ah); 0250 C = y/Min + a*Ah*Inv*b; 0251 Ch = x/Max + ah*Inv*A*bh; 0252 d1 = Qr(i,i+1)*(1-Qrh(i+1,i)); 0253 D = d1+a*eh+a*Ah*Inv*A*eh; 0254 d1h = Qrh(j,j-1)*(1-Qr(j-1,j)); 0255 Dh = d1h+ah*Inv*e; 0256 d0 = sum(Qr(i,i+1:j-1)); 0257 E = 1-d0-y/Min+a*Ah*Inv*e; 0258 d0h = sum(Qrh(j,i+1:j-1)); 0259 Eh = 1-d0h-x/Max+ah*Inv*A*eh; 0260 r1 = Qr(i,i+1)*Qrh(i+1,i); 0261 R = r1+a*bh+a*Ah*Inv*A*bh; 0262 r1h = Qrh(j,j-1)*Qr(j-1,j); 0263 Rh = r1h+ah*Inv*b; 0264 0265 % For Part 3 of cycle 0266 0267 A3 = Q(I,J); 0268 A3h = Qh(J,I); 0269 a3 = Q(i,J); 0270 a3h = Qh(j,I); 0271 e3 = 1 - sum(Q(I,i+2:j-1),2); 0272 e3h = 1 - sum(Qh(J,i+1:j-2),2); 0273 0274 Inv3 = inv(Eye-A3*A3h); 0275 d3 = sum(Q(i,i+1:j-1)); 0276 E3 = 1-d3 + a3*A3h*Inv3*e3; 0277 d3h = sum(Qh(j,i+1:j-1)); 0278 E3h = 1-d3h + a3h*Inv3*A3*e3h; 0279 0280 end 0281 0282 if ~(i == 1 & j == n) 0283 0284 % Standing 0285 if j == n 0286 pS1 = 0; pS2=0; pS3=0; 0287 else 0288 0289 % Part 1 of cycle 0290 ES = E + C*Dh/(1-Rh); 0291 RS = R + C*Ch/(1-Rh); 0292 pS1 = ES/(1-RS); 0293 0294 % Part 2 of cycle 0295 pS2 = Ch/(1-Rh); 0296 0297 % Part 3 of cycle 0298 pS3 = E3h; 0299 end 0300 0301 % Hanging 0302 if i == 1 0303 pH1=0; pH2=0; pH3=0; 0304 else 0305 0306 % Part 1 of cycle 0307 EH = Eh + Ch*D/(1-R); 0308 RH = Rh + Ch*C/(1-R); 0309 pH1 = EH/(1-RH); 0310 0311 % Part 2 of cycle 0312 pH2 = C/(1-R); 0313 0314 % Part 3 of cycle 0315 pH3 = E3; 0316 end 0317 0318 % Check 0319 0320 end 0321 0322 if Min == 0 | Max == 0 0323 pS = 0; 0324 pH = 0; 0325 elseif ~(i == 1 & j == n) 0326 pS = Max*pS1*pS2*pS3; 0327 pH = Min*pH1*pH2*pH3; 0328 else % i == 1 & j == n 0329 pS = Min*E3; % = Max*E3h; 0330 pH = 0; 0331 % Old code 0332 % pS = 1/2*Min*E3; 0333 % pH = 1/2*Max*E3h; 0334 end 0335 0336 F_rfc(i,j) = pS; % Standing 0337 F_rfc(j,i) = pH; % Hanging 0338 0339 end 0340 end 0341 %fprintf(1,'\n'); 0342 0343 % Check negative elements 0344 0345 [I,J] = find(F_rfc<0); 0346 if ~isempty(I) 0347 warning(['Negative elements in F_rfc. Setting to zero!']); 0348 for k = 1: length(I) 0349 F_rfc(I(k),J(k)) = 0; 0350 end 0351 end 0352 0353 % Multiply with the intensity of local minima for the swithching process 0354 0355 F_rfc = c_m*F_rfc; 0356 0357 T(6,:)=clock; 0358 0359 0360
Comments or corrections to the WAFO group