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

## CROSS-REFERENCE INFORMATION

This function calls:
 clock Current date and time as date vector. error Display message and abort function. inv Matrix inverse. tril Extract lower triangular part. triu Extract upper triangular part.
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 %
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