Home > wafo > markov > arfm2mctp.m

arfm2mctp

PURPOSE ^

Calculates the markov matrix given an asymmetric rainflow matrix.

SYNOPSIS ^

[F,T] = arfm2mctp(Frfc)

DESCRIPTION ^

 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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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

Mathematical Statistics
Centre for Mathematical Sciences
Lund University with Lund Institute of Technology

Comments or corrections to the WAFO group


Generated on Thu 06-Oct-2005 02:21:16 for WAFO by m2html © 2003