ARFM2MCTP Calculates the markov matrix given an asymmetric rainflow matrix. CALL: F = arfm2mctp(Frfc); F = Markov matrix (from-to-matrix) [n,n] Frfc = Rainflow Matrix [n,n] Example: param = [-1 1 32]; u = levels(param); F = mktestmat(param,[-0.2 0.2],0.15,2); F = F/sum(sum(F)); Farfc = mctp2arfm({F []}); F1 = arfm2mctp(Farfc); cmatplot(u,u,{F+F' F1},3); sum(sum(abs(F1-(F+F')))) % should be zero See also rfm2mctp, mctp2arfm, smctp2arfm, cmatplot
Current date and time as date vector. | |
Display message and abort function. | |
Matrix inverse. | |
Extract lower triangular part. | |
Extract upper triangular part. |
0001 function [F,T] = arfm2mctp(Frfc) 0002 %ARFM2MCTP Calculates the markov matrix given an asymmetric rainflow matrix. 0003 % 0004 % CALL: F = arfm2mctp(Frfc); 0005 % 0006 % F = Markov matrix (from-to-matrix) [n,n] 0007 % Frfc = Rainflow Matrix [n,n] 0008 % 0009 % Example: 0010 % param = [-1 1 32]; u = levels(param); 0011 % F = mktestmat(param,[-0.2 0.2],0.15,2); 0012 % F = F/sum(sum(F)); 0013 % Farfc = mctp2arfm({F []}); 0014 % F1 = arfm2mctp(Farfc); 0015 % cmatplot(u,u,{F+F' F1},3); 0016 % sum(sum(abs(F1-(F+F')))) % should be zero 0017 % 0018 % See also rfm2mctp, mctp2arfm, smctp2arfm, cmatplot 0019 0020 % References: 0021 % 0022 % P. Johannesson (1999): 0023 % Rainflow Analysis of Switching Markov Loads. 0024 % PhD thesis, Mathematical Statistics, Centre for Mathematical Sciences, 0025 % Lund Institute of Technology. 0026 0027 % Tested on Matlab 5.3 0028 % 0029 % History: 0030 % Revised by PJ 09-Apr-2001 0031 % updated for WAFO 0032 % Created by PJ (Pär Johannesson) 1998 0033 % Copyright (c) 1997-1998 by Pär Johannesson 0034 % Toolbox: Rainflow Cycles for Switching Processes V.1.1, 22-Jan-1998 0035 0036 % Recursive formulation a'la Igor 0037 % 0038 % This program used the formulation where the probabilities 0039 % of the events are calculated using "elementary" events for 0040 % the MCTP. 0041 % 0042 % Standing 0043 % pS = Max*pS1*pS2*pS3; 0044 % F_rfc(i,j) = pS; 0045 % Hanging 0046 % pH = Min*pH1*pH2*pH3; 0047 % F_rfc(j,i) = pH; 0048 % 0049 % The cond. prob. pS1, pS2, pS3, pH1, pH2, pH3 are calculated using 0050 % the elementary cond. prob. C, E, R, D, E3, Ch, Eh, Rh, Dh, E3h. 0051 0052 ni = nargin; 0053 no = nargout; 0054 error(nargchk(1,1,ni)); 0055 0056 T(1,:)=clock; 0057 0058 N = sum(sum(Frfc)); 0059 Frfc = Frfc/N; 0060 0061 n = length(Frfc); % Number of levels 0062 0063 T(7,:)=clock; 0064 0065 % Transition matrices for MC 0066 0067 Q = zeros(n,n); 0068 Qh = zeros(n,n); 0069 0070 % Transition matrices for time-reversed MC 0071 0072 Qr = zeros(n,n); 0073 Qrh = zeros(n,n); 0074 0075 % Probability of minimum and of maximun 0076 0077 MIN = sum(triu(Frfc)') + sum(tril(Frfc)); 0078 MAX = sum(triu(Frfc)) + sum(tril(Frfc)'); 0079 0080 % Calculate rainflow matrix 0081 0082 F = zeros(n,n); 0083 EYE = eye(n,n); 0084 0085 %fprintf(1,'Calculating row '); 0086 for k=1:n-1 % k = subdiagonal 0087 % fprintf(1,'-%1d',i); 0088 0089 for i=1:n-k % i = minimum 0090 0091 j = i+k; % maximum; 0092 0093 pS = Frfc(i,j); % Standing cycle 0094 pH = Frfc(j,i); % Hanging cycle 0095 0096 Min = MIN(i); 0097 Max = MAX(j); 0098 0099 % fprintf(1,'Min=%f, Max=%f\n',Min,Max); 0100 0101 0102 if j-i == 2 % Second subdiagonal 0103 0104 % For Part 1 & 2 of cycle 0105 0106 %C = y/Min; 0107 c0 = 0; 0108 c1 = 1/Min; 0109 %Ch = x/Max; 0110 c0h = 0; 0111 c1h = 1/Max; 0112 d1 = Qr(i,i+1)*(1-Qrh(i+1,i)); 0113 D = d1; 0114 d1h = Qrh(j,j-1)*(1-Qr(j-1,j)); 0115 Dh = d1h; 0116 d0 = sum(Qr(i,i+1:j-1)); 0117 %E = 1-d0-y/Min; 0118 e0 = 1-d0; 0119 e1 = -1/Min; 0120 d0h = sum(Qrh(j,i+1:j-1)); 0121 %Eh = 1-d0h-x/Max; 0122 e0h = 1-d0h; 0123 e1h = -1/Max; 0124 r1 = Qr(i,i+1)*Qrh(i+1,i); 0125 R = r1; 0126 r1h = Qrh(j,j-1)*Qr(j-1,j); 0127 Rh = r1h; 0128 0129 % For Part 3 of cycle 0130 0131 d3h = sum(Qh(j,i+1:j-1)); 0132 E3h = 1-d3h; 0133 d3 = sum(Q(i,i+1:j-1)); 0134 E3 = 1-d3; 0135 0136 % Define coeficients for equation system 0137 a0 = -pS+2*pS*Rh-pS*Rh^2+pS*R-2*pS*Rh*R+pS*Rh^2*R; 0138 a1 = -E3h*Max*c1h*e0*Rh+E3h*Max*c1h*e0; 0139 a3 = -E3h*Max*c1h*e1*Rh+E3h*Max*c1h*Dh*c1+E3h*Max*c1h*e1+pS*c1h*c1-pS*c1h*c1*Rh; 0140 0141 b0 = -pH+2*pH*R+pH*Rh-2*pH*Rh*R-pH*R^2+pH*Rh*R^2; 0142 b2 = -Min*E3*e0h*R*c1+Min*E3*e0h*c1; 0143 b3 = Min*E3*e1h*c1+Min*E3*D*c1h*c1-pH*c1h*c1*R-Min*E3*e1h*R*c1+pH*c1h*c1; 0144 0145 C2 = a3*b2; 0146 C1 = (-a0*b3+a1*b2+a3*b0); 0147 C0 = a1*b0; 0148 % Solve: C2*z^2 + C1*z + C0 = 0 0149 z1 = -C1/2/C2 + sqrt((C1/2/C2)^2-C0/C2); 0150 z2 = -C1/2/C2 - sqrt((C1/2/C2)^2-C0/C2); 0151 0152 % Solution 1 0153 x1 = -(b0+b2*z1)/(b3*z1); 0154 y1 = z1; 0155 % Solution 2 0156 x2 = -(b0+b2*z2)/(b3*z2); 0157 y2 = z2; 0158 0159 x = x2; 0160 y = y2; 0161 0162 % fprintf(1,'2nd: i=%d, j=%d: x1=%f, y1=%f, x2=%f, y2=%f\n',i,j,x1,y1,x2,y2); 0163 0164 % Test Standing cycle: assume x=y 0165 0166 C0 = a0; C1 = a1; C2 = a3; 0167 z1S = -C1/2/C2 + sqrt((C1/2/C2)^2-C0/C2); 0168 z2S = -C1/2/C2 - sqrt((C1/2/C2)^2-C0/C2); 0169 0170 % Test Hanging cycle: assume x=y 0171 0172 C0 = b0; C1 = b2; C2 = b3; 0173 z1H = -C1/2/C2 + sqrt((C1/2/C2)^2-C0/C2); 0174 z2H = -C1/2/C2 - sqrt((C1/2/C2)^2-C0/C2); 0175 0176 % fprintf(1,'2nd: i=%d, j=%d: z1S=%f,: z2S=%f, z1H=%f, z2H=%f\n',i,j,z1S,z2S,z1H,z2H); 0177 0178 else 0179 0180 Eye = EYE(1:j-i-2,1:j-i-2); 0181 0182 % For Part 1 & 2 of cycle 0183 0184 I = i+1:j-2; 0185 J = i+2:j-1; 0186 A = Qr(I,J); 0187 Ah = Qrh(J,I); 0188 a = Qr(i,J); 0189 ah = Qrh(j,I); 0190 b = Qr(I,j); 0191 bh = Qrh(J,i); 0192 0193 e = 1 - sum(Qr(I,i+2:j),2); 0194 eh = 1 - sum(Qrh(J,i:j-2),2); 0195 0196 Inv = inv(Eye-A*Ah); 0197 %C = y/Min + a*Ah*Inv*b; 0198 c0 = a*Ah*Inv*b; 0199 c1 = 1/Min; 0200 %Ch = x/Max + ah*Inv*A*bh; 0201 c0h = ah*Inv*A*bh; 0202 c1h = 1/Max; 0203 d1 = Qr(i,i+1)*(1-Qrh(i+1,i)); 0204 D = d1+a*eh+a*Ah*Inv*A*eh; 0205 d1h = Qrh(j,j-1)*(1-Qr(j-1,j)); 0206 Dh = d1h+ah*Inv*e; 0207 d0 = sum(Qr(i,i+1:j-1)); 0208 %E = 1-d0-y/Min+a*Ah*Inv*e; 0209 e0 = 1-d0+a*Ah*Inv*e; 0210 e1 = -1/Min; 0211 d0h = sum(Qrh(j,i+1:j-1)); 0212 %Eh = 1-d0h-x/Max+ah*Inv*A*eh; 0213 e0h = 1-d0h+ah*Inv*A*eh; 0214 e1h = -1/Max; 0215 r1 = Qr(i,i+1)*Qrh(i+1,i); 0216 R = r1+a*bh+a*Ah*Inv*A*bh; 0217 r1h = Qrh(j,j-1)*Qr(j-1,j); 0218 Rh = r1h+ah*Inv*b; 0219 0220 % For Part 3 of cycle 0221 0222 A3 = Q(I,J); 0223 A3h = Qh(J,I); 0224 Inv3 = inv(Eye-A3*A3h); 0225 0226 % For Standing cycle 0227 d3h = sum(Qh(j,i+1:j-1)); 0228 c3h = Qh(j,I); 0229 e3h = 1 - sum(Qh(J,i+1:j-2),2); 0230 E3h = 1-d3h + c3h*Inv3*A3*e3h; 0231 0232 % For Hanging cycle 0233 d3 = sum(Q(i,i+1:j-1)); 0234 c3 = Q(i,J); 0235 e3 = 1 - sum(Q(I,i+2:j-1),2); 0236 E3 = 1-d3 + c3*A3h*Inv3*e3; 0237 0238 end 0239 0240 if j-i == 1 % First subdiagonal 0241 0242 if i == 1 0243 x = Max; 0244 y = Max; 0245 elseif j == n 0246 x = Min; 0247 y = Min; 0248 else 0249 if pS == 0 0250 x = 0; 0251 y = pH; 0252 elseif pH == 0 0253 x = pS; 0254 y = 0; 0255 else 0256 x = Min*pS/(Min-pH); 0257 y = Max*pH/(Max-pS); 0258 end 0259 end 0260 0261 elseif j-i >= 2 0262 if i == 1 0263 x = Max*(1-sum(Qh(j,2:j-1))); 0264 y = Max*(1-sum(Qrh(j,2:j-1))); 0265 elseif j == n 0266 x = Min*(1-sum(Q(i,i+1:n-1))); 0267 y = Min*(1-sum(Qr(i,i+1:n-1))); 0268 else 0269 if pS == 0 0270 x = 0; 0271 y = pH; 0272 elseif pH == 0 0273 x = pS; 0274 y = 0; 0275 else 0276 % Define coeficients for equation system 0277 a0 = pS*c0h*c0+pS*Rh^2*R-2*pS*Rh*R-E3h*Max*c0h*e0*Rh+E3h*Max*c0h*e0+2*pS*Rh+pS*R-pS*c0h*c0*Rh-pS-pS*Rh^2+E3h*Max*c0h*Dh*c0; 0278 a1 = pS*c1h*c0+E3h*Max*c1h*Dh*c0-E3h*Max*c1h*e0*Rh-pS*c1h*c0*Rh+E3h*Max*c1h*e0; 0279 a2 = pS*c0h*c1+E3h*Max*c0h*e1-pS*c0h*c1*Rh+E3h*Max*c0h*Dh*c1-E3h*Max*c0h*e1*Rh; 0280 a3 = -E3h*Max*c1h*e1*Rh+E3h*Max*c1h*Dh*c1+E3h*Max*c1h*e1+pS*c1h*c1-pS*c1h*c1*Rh; 0281 0282 b0 = pH*c0h*c0+pH*Rh*R^2-pH+pH*Rh-2*pH*Rh*R-pH*c0h*c0*R+Min*E3*e0h*c0-Min*E3*e0h*R*c0+Min*E3*D*c0h*c0+2*pH*R-pH*R^2; 0283 b1 = Min*E3*D*c1h*c0+Min*E3*e1h*c0+pH*c1h*c0-Min*E3*e1h*R*c0-pH*c1h*c0*R; 0284 b2 = -pH*c0h*c1*R-Min*E3*e0h*R*c1+Min*E3*D*c0h*c1+Min*E3*e0h*c1+pH*c0h*c1; 0285 b3 = Min*E3*e1h*c1+Min*E3*D*c1h*c1-pH*c1h*c1*R-Min*E3*e1h*R*c1+pH*c1h*c1; 0286 0287 C2 = a2*b3-a3*b2; 0288 C1 = a0*b3-a1*b2+a2*b1-a3*b0; 0289 C0 = a0*b1-a1*b0; 0290 %fprintf(1,'i=%d, j=%d, C0/C2=%f,C1/C2=%f,C2=%f\n',i,j,C0/C2,C1/C2,C2); 0291 % Solve: C2*z^2 + C1*z + C0 = 0 0292 z1 = -C1/2/C2 + sqrt((C1/2/C2)^2-C0/C2); 0293 z2 = -C1/2/C2 - sqrt((C1/2/C2)^2-C0/C2); 0294 0295 % Solution 1 0296 x1 = -(b0+b2*z1)/(b1+b3*z1); 0297 y1 = z1; 0298 % Solution 2 0299 x2 = -(b0+b2*z2)/(b1+b3*z2); 0300 y2 = z2; 0301 0302 x = x2; 0303 y = y2; 0304 0305 % fprintf(1,'End: i=%d, j=%d: x1=%f, y1=%f, x2=%f, y2=%f\n',i,j,x1,y1,x2,y2); 0306 end 0307 end 0308 end 0309 0310 % fprintf(1,'i=%d, j=%d: x=%f, y=%f\n',i,j,x,y); 0311 0312 % min-max 0313 F(i,j) = x; 0314 0315 % max-min 0316 F(j,i) = y; 0317 0318 % Fill the transitions matrices 0319 Q(i,j) = x/Min; 0320 Qh(j,i) = y/Max; 0321 Qr(i,j) = y/Min; 0322 Qrh(j,i) = x/Max; 0323 0324 end 0325 end 0326 %fprintf(1,'\n'); 0327 0328 0329 T(8,:)=clock; 0330 0331 0332 0333 0334
Comments or corrections to the WAFO group