Home > wafo > markov > mctp2arfm.m

# mctp2arfm

## PURPOSE

Calculates asymmetric rainflow matrix for a MCTP.

## SYNOPSIS

[F_rfc,FF,FFr,T] = mctp2arfm(F,c_m)

## DESCRIPTION

``` 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

## CROSS-REFERENCE INFORMATION

This function calls:
 mat2tmat Converts a matrix to a transition matrix. mc2stat Calculates the stationary distribution for a Markov chain. clock Current date and time as date vector. error Display message and abort function. inv Matrix inverse. warning Display warning message; disable or enable warning messages.
This function is called by:
 test_markov Quick test of the routines in module 'markov'

## SOURCE CODE

```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```

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