Home > wafo > markov > smctp2arfm.m

smctp2arfm

PURPOSE ^

Calculates the asymmetric rainflow matrix for a SMCTP.

SYNOPSIS ^

[F_rfc,F_rfc_z,T] = smctp2arfm(P,F,c_m,SideInfo)

DESCRIPTION ^

 SMCTP2ARFM  Calculates the asymmetric rainflow matrix for a SMCTP.
 
  CALL:  [Frfc]        = smctp2arfm(P,F);
         [Frfc,Frfc_z] = smctp2arfm(P,F,c_m,SideInfo);
 
  Frfc    = Rainflow Matrix                            [n,n]
  Frfc_z  = Rainflow Matrix, with side info z    {r,r1}[n,n]
 
  P       = Transition matrix for regime process.      [r,r]
  F       = Cell array of min-Max and Max-min matrices {r,2}
  F{i,1}  = min-Max matrix, process i                  [n,n]
  F{i,2}  = Max-min matrix, process i                  [n,n]
  c_m     = Intensity of local minima, switching proc. [1,1]
            (Default: 1)
  SideInfo = Which type of side information
           0: No side info
           1: Mark min & max, r1=r
           2: Mark when counted, r1=1
 
  Calculates the asymmetric rainflow matrix for a switching process
  with a Markov chain of turning points within each regime.
  If a matrix F{i,2}=[], then the process will be assumed to be 
  time-reversible.
 
  Example: (Two processes as in Example 4.1 in PhD thesis)
    P = [0.9 0.1; 0.05 0.95];
    param = [-1 1 32]; u = levels(param);
    F1 = mktestmat(param,[-0.4 -0.3],0.15,1);
    F2 = mktestmat(param,[0.3 0.4],0.15,1);
    [Frfc,Frfc_z] = smctp2arfm(P,{F1 F1'; F2 F2'},1,1);
    figure(1),cmatplot(u,u,Frfc_z,3)
    Frfc1 = Frfc_z{1,1}+Frfc_z{1,2}+Frfc_z{2,1}+Frfc_z{2,2};
    figure(2),cmatplot(u,u,{Frfc Frfc1},3) % Shall be identical
 
  See also  mctp2arfm, smctp2rfm, dtp2arfm_sid

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [F_rfc,F_rfc_z,T] = smctp2arfm(P,F,c_m,SideInfo)
0002 %SMCTP2ARFM  Calculates the asymmetric rainflow matrix for a SMCTP.
0003 %
0004 % CALL:  [Frfc]        = smctp2arfm(P,F);
0005 %        [Frfc,Frfc_z] = smctp2arfm(P,F,c_m,SideInfo);
0006 %
0007 % Frfc    = Rainflow Matrix                            [n,n]
0008 % Frfc_z  = Rainflow Matrix, with side info z    {r,r1}[n,n]
0009 %
0010 % P       = Transition matrix for regime process.      [r,r]
0011 % F       = Cell array of min-Max and Max-min matrices {r,2}
0012 % F{i,1}  = min-Max matrix, process i                  [n,n]
0013 % F{i,2}  = Max-min matrix, process i                  [n,n]
0014 % c_m     = Intensity of local minima, switching proc. [1,1]
0015 %           (Default: 1)
0016 % SideInfo = Which type of side information
0017 %          0: No side info
0018 %          1: Mark min & max, r1=r
0019 %          2: Mark when counted, r1=1
0020 %
0021 % Calculates the asymmetric rainflow matrix for a switching process
0022 % with a Markov chain of turning points within each regime.
0023 % If a matrix F{i,2}=[], then the process will be assumed to be 
0024 % time-reversible.
0025 %
0026 % Example: (Two processes as in Example 4.1 in PhD thesis)
0027 %   P = [0.9 0.1; 0.05 0.95];
0028 %   param = [-1 1 32]; u = levels(param);
0029 %   F1 = mktestmat(param,[-0.4 -0.3],0.15,1);
0030 %   F2 = mktestmat(param,[0.3 0.4],0.15,1);
0031 %   [Frfc,Frfc_z] = smctp2arfm(P,{F1 F1'; F2 F2'},1,1);
0032 %   figure(1),cmatplot(u,u,Frfc_z,3)
0033 %   Frfc1 = Frfc_z{1,1}+Frfc_z{1,2}+Frfc_z{2,1}+Frfc_z{2,2};
0034 %   figure(2),cmatplot(u,u,{Frfc Frfc1},3) % Shall be identical
0035 %
0036 % See also  mctp2arfm, smctp2rfm, dtp2arfm_sid
0037 
0038 % References:
0039 %  
0040 %  P. Johannesson (1999):
0041 %  Rainflow Analysis of Switching Markov Loads.
0042 %  PhD thesis, Mathematical Statistics, Centre for Mathematical Sciences,
0043 %  Lund Institute of Technology.
0044   
0045 % Tested  on Matlab  5.3
0046 %
0047 % History:
0048 % Revised by PJ  09-Apr-2001
0049 %   updated for WAFO
0050 % Created by PJ (Pär Johannesson) 1998
0051 % Copyright (c) 1997-1998 by Pär Johannesson
0052 % Toolbox: Rainflow Cycles for Switching Processes V.1.1, 22-Jan-1998
0053 
0054 
0055 % This program used the formulation where the probabilities
0056 % of the events are calculated using "elementary" events for
0057 % the MCTP.
0058 %
0059 % Standing
0060 %    pS = Max*pS1*pS2*pS3;
0061 %    F_rfc(i,j) = pS;
0062 % Hanging
0063 %    pH = Min*pH1*pH2*pH3;
0064 %    F_rfc(j,i) = pH;
0065 %
0066 % The cond. prob. pS1, pS2, pS3, pH1, pH2, pH3 are calculated using
0067 % the elementary cond. prob. C, E, R, D, E3, Ch, Eh, Rh, Dh, E3h. 
0068 
0069 T(1,:)=clock;
0070 
0071 % Check input arguments
0072 
0073 ni = nargin;
0074 no = nargout;
0075 error(nargchk(2,4,ni));
0076 
0077 if ni < 3
0078   c_m=[];
0079 end
0080 
0081 if ni < 4
0082   SideInfo = [];
0083 end
0084 
0085 if isempty(c_m)
0086   c_m=1;
0087 end
0088 if isempty(SideInfo)
0089   SideInfo=0;
0090 end
0091 
0092 % Define 
0093 
0094 Zstr = '123456789'; Zstr=Zstr';
0095 
0096 r = length(P);   % Number of regime states
0097 n = length(F{1,1});  % Number of levels
0098 
0099 % Check that the rowsums of P are equal to 1
0100 
0101 P = mat2tmat(P);
0102 
0103 T(2,:)=clock;
0104 
0105 % Normalize the rowsums of F{1,1},...,F{r,1} to 1
0106 %  ==>  Q{1,1},...,Q{r,1}
0107 
0108 for i = 1:r
0109   QQ{i,1} = F{i,1};
0110   QQ{i,1} = mat2tmat(QQ{i,1},1);
0111 end
0112 
0113 T(3,:)=clock;
0114 
0115 % Normalize the rowsums of F{1,2},...,F{r,2} to 1
0116 %  ==>  Q{1,2},...,Q{r,2}
0117 
0118 for i = 1:r
0119   
0120   if isempty(F{i,2})        % Time-reversible
0121     QQ{i,2} = F{i,1};
0122     QQ{i,2} = QQ{i,2}';  
0123   else                   % Fhi is given
0124     QQ{i,2} = F{i,2}; 
0125   end
0126     
0127   QQ{i,2} = mat2tmat(QQ{i,2},-1);
0128 
0129 end
0130 
0131 T(4,:)=clock;
0132 
0133 % Make the transition matrix Q for the joint min-Max process
0134 
0135 Q = zeros(n*r,n*r);
0136 I = 0:r:(n-1)*r;
0137 for z=1:r
0138   Q0 = kron(QQ{z,1},P);
0139   Q(I+z,:) = Q0(I+z,:);
0140 end
0141 
0142 T(5,:)=clock;
0143 
0144 % Make the transition matrix Qh for the joint Max-min process
0145 
0146 Qh = zeros(n*r,n*r);
0147 I = 0:r:(n-1)*r;
0148 for z=1:r
0149   Q0 = kron(QQ{z,2},P);
0150   Qh(I+z,:) = Q0(I+z,:);
0151 end
0152 
0153 T(6,:)=clock;
0154 
0155 % Stationary distribution (=ro) of local minima with transition matrix
0156 % Qt = Q*Qh = "Transition matrix for min-to-min"
0157 
0158 Qt = Q*Qh;
0159 ro = mc2stat(Qt(1:r*(n-1),1:r*(n-1)));  % Stationary distr., row vector  
0160 ro = [ro zeros(1,r)];  % Minimum can't reach the highest level
0161 
0162 % Stationary distribution (=roh) of local maxima with transition matrix
0163 % Qt = Qh*Q = "Transition matrix for max-to-max"
0164 
0165 Qth = Qh*Q;
0166 roh = mc2stat(Qth(r+1:r*(n),r+1:r*(n)));  % Stationary distr., row vector  
0167 roh = [zeros(1,r) roh];  % Maximum can't reach the highest level
0168 
0169 T(7,:)=clock;
0170 
0171 % Make the frequency matrix FF for the joint min-Max and Max-min
0172 % distribution
0173 
0174 FF = Q.*(ro'*ones(1,n*r)) + Qh.*(roh'*ones(1,n*r));
0175 
0176 % Create Transition matrices for time-reversed MC
0177 
0178 % Backward min-to-max
0179 I1 = find(ro>0); I2 = find(ro<=0);
0180 ro_inv = ro; ro_inv(I1) = 1./ro(I1); ro_inv(I2) = zeros(1,length(I2));
0181 Qr = Qh' .* (ro_inv'*roh);
0182 
0183 % Backward max-to-min
0184 I1 = find(roh>0); I2 = find(roh<=0);
0185 roh_inv = roh; roh_inv(I1) = 1./roh(I1); roh_inv(I2) = zeros(1,length(I2));
0186 Qrh = Q' .* (roh_inv'*ro);
0187 
0188 % Make the frequency matrix FF for the joint min-Max and Max-min
0189 % distribution
0190 
0191 FF1 = Qr.*(ro'*ones(1,n*r)) + Qrh.*(roh'*ones(1,n*r));
0192 
0193 T(8,:)=clock;
0194 
0195 % Initiation of matrices
0196 
0197 F_rfc = zeros(n,n);
0198 EYE = eye(n*r,n*r); Eye1  = eye(r,r);
0199 Zero1 = zeros(r,1); Zero2 = zeros(r,r);
0200 One1 = ones(r,1);   One2  = ones(r,r);
0201 
0202 if SideInfo == 1  % Rainflow matrix with side information
0203   [F_rfc_z{1:r,1:r}] = deal(zeros(n,n));
0204 end % Rainflow matrix with side information
0205     
0206 % Calculate rainflow matrix 
0207 
0208 
0209 %fprintf(1,'Calculating row ');
0210 for i=1:n-1
0211 %  fprintf(1,'-%1d',i);
0212 
0213   for j=i+1:n
0214 
0215     I = r*(i-1)+1:r*i;
0216     J = r*(j-1)+1:r*j;
0217     I1 = r*(i+1-1)+1:r*(i+1);
0218     J1 = r*(j-1-1)+1:r*(j-1);
0219     I1J1 = r*(i+1-1)+1:r*(j-1);
0220 
0221     x = FF(I,J);
0222     y = FF(J,I);
0223 
0224     Min = sum(FF(I,:),2)';
0225     Max = sum(FF(J,:),2)';
0226 
0227     Ro = ro(I);    % Probability of "min=i"
0228     Roh = roh(J);  % Probability of "max=j"
0229 %    fprintf(1,'Min=%f, Max=%f, Ro=%f, Roh=%f\n',Min,Max,Ro,Roh);
0230 
0231     Min = Ro; Max = Roh; % Just to be sure
0232     
0233     Min = Min'*One1';
0234     Max = Max'*One1';
0235     
0236     % For all subdiagonals
0237 
0238     C = y'./Min; Ch = x'./Max;
0239     E = 1-sum(y'./Min,2); Eh = 1-sum(x'./Max,2);
0240 
0241     
0242     if j-i >= 2  % Second and higher subdiagonals
0243 
0244       % For Part 1 & 2 of cycle
0245 
0246       d1  = Qr(I,I1)*(1-sum(Qrh(I1,I),2));
0247       D   = d1;
0248       d1h = Qrh(J,J1)*(1-sum(Qr(J1,J),2));
0249       Dh  = d1h;
0250       d0  = sum(Qr(I,I1J1),2); 
0251       E   = E-d0; 
0252       d0h = sum(Qrh(J,I1J1),2); % ???
0253       Eh  = Eh-d0h;
0254       r1  = Qr(I,I1)*Qrh(I1,I);
0255       R   = r1;
0256       r1h = Qrh(J,J1)*Qr(J1,J);
0257       Rh  = r1h;
0258 
0259       % For Part 3 of cycle
0260 
0261       d3  = sum(Q(I,I1J1),2);
0262       E3  = 1-d3;
0263       d3h = sum(Qh(J,I1J1),2);
0264       E3h = 1-d3h;
0265 
0266     else % First subdiagonal
0267       
0268       D  = Zero1; Dh  = Zero1;
0269       R  = Zero2; Rh  = Zero2;
0270       E3 = One1;  E3h = One1;
0271       
0272     end
0273     
0274       
0275     if j-i >= 3  % Third and higher subdiagonals
0276 
0277       Eye = EYE(1:r*(j-i-2),1:r*(j-i-2));
0278 
0279       % For Part 1 & 2 of cycle
0280 
0281       II  = r*(i+1-1)+1:r*(j-2);  % i+1:j-2
0282       JJ  = r*(i+2-1)+1:r*(j-1);  % i+2:j-1;
0283       
0284       A  = Qr(II,JJ);
0285       Ah = Qrh(JJ,II);
0286       Inv = inv(Eye-A*Ah);
0287       
0288       a  = Qr(I,JJ);
0289       ah = Qrh(J,II);
0290       b  = Qr(II,J);
0291       bh = Qrh(JJ,I);
0292 
0293       e  = 1 - sum(Qr(II,r*(i+2-1)+1:r*(j)),2);   % i+2:j
0294       eh = 1 - sum(Qrh(JJ,r*(i-1)+1:r*(j-2)),2);  % i:j-2
0295 
0296       C   = C + a*Ah*Inv*b;
0297       Ch  = Ch + ah*Inv*A*bh;
0298       D   = D + a*eh+a*Ah*Inv*A*eh;
0299       Dh  = Dh + ah*Inv*e;
0300       E   = E + a*Ah*Inv*e;
0301       Eh  = Eh + ah*Inv*A*eh;
0302       R   = R + a*bh+a*Ah*Inv*A*bh;
0303       Rh  = Rh + ah*Inv*b;
0304 
0305       % For Part 3 of cycle
0306 
0307       A3  = Q(II,JJ);
0308       A3h = Qh(JJ,II);
0309       Inv3 = inv(Eye-A3*A3h);
0310       c3  = Q(I,JJ);
0311       c3h = Qh(J,II);
0312       e3  = 1 - sum(Q(II,r*(i+2-1)+1:r*(j-1)),2);   % i+2:j-1
0313       e3h  = 1 - sum(Qh(JJ,r*(i+1-1)+1:r*(j-2)),2); % i+1:j-2
0314 
0315       E3  = E3  + c3*A3h*Inv3*e3;
0316       E3h = E3h + c3h*Inv3*A3*e3h;
0317 
0318     end
0319 
0320     if ~(i == 1 & j == n)
0321 
0322       % Standing
0323       if j == n
0324         pS1 = Zero1; pS2=Zero2; pS3=Zero1;
0325       else
0326 
0327         % Part 1 of cycle
0328     AS = (Eye1-Rh)\Dh;
0329         ES = E + C*AS;
0330     BS = (Eye1-Rh)\Ch;
0331         RS = R + C*BS;
0332         pS1 = (Eye1-RS)\ES;
0333 
0334         % Part 2 of cycle
0335         pS2 = (Eye1-Rh)\Ch;
0336 
0337         % Part 3 of cycle
0338         pS3 = E3h;
0339     
0340       end
0341 
0342       % Hanging
0343       if i == 1
0344         pH1 = Zero1; pH2=Zero2; pH3=Zero1;
0345       else
0346 
0347         % Part 1 of cycle
0348     AH = (Eye1-R)\D;
0349         EH = Eh + Ch*AH;
0350     BH = (Eye1-R)\C;
0351         RH = Rh + Ch*BH;
0352         pH1 = (Eye1-RH)\EH;
0353 
0354         % Part 2 of cycle
0355         pH2 = (Eye1-R)\C;
0356 
0357         % Part 3 of cycle
0358         pH3 = E3;
0359       end
0360 
0361     else % i == 1 & j == n
0362       
0363       % Standing
0364       pS1 = One1; 
0365       pS2 = Eye1;
0366       pS3 = E3h*sum(Roh)/sum(Ro+Roh);
0367       
0368       % Hanging
0369       pH1 = One1; 
0370       pH2 = Eye1;
0371       pH3 = E3*sum(Ro)/sum(Ro+Roh);
0372       
0373     end
0374 
0375     % Standing RFC
0376     F_rfc(i,j) = Roh*diag(pS3)*pS2*pS1;
0377     
0378     % Hanging RFC
0379     F_rfc(j,i) = Ro*diag(pH3)*pH2*pH1;
0380     
0381     if SideInfo == 1  % Rainflow matrix with side information
0382       if (i==1) & (j==n)
0383     
0384         b3  = Q(II,J);
0385         EE  = Q(I,J) + c3*A3h*Inv3*b3;
0386         for z=1:r
0387         for w=1:r
0388       
0389             % Standing RFC
0390         F_rfc_z{z,w}(i,j) = 1/2*Ro(z)*EE(z,w);
0391       
0392             % Hanging RFC
0393         F_rfc_z{z,w}(j,i) = 1/2*Ro(z)*EE(z,w);
0394       end
0395     end
0396     
0397       else
0398     
0399         for z=1:r
0400         for w=1:r
0401       
0402             % Standing RFC
0403         F_rfc_z{z,w}(i,j) = Roh(w)*pS3(w)*pS2(w,z)*pS1(z);
0404       
0405             % Hanging RFC
0406         F_rfc_z{z,w}(j,i) = Ro(w)*pH3(w)*pH2(w,z)*pH1(z);
0407       
0408       end
0409     end
0410     
0411       end
0412     elseif SideInfo ~= 0
0413       error(['SideInfo = ' num2str(SideInfo) ' not implemented']);
0414     end % Rainflow matrix with side information
0415     
0416   end 
0417 end
0418 %fprintf(1,'\n');
0419 
0420 % Multiply with the intensity of local minima for the swithching process
0421 
0422 F_rfc = c_m*F_rfc;
0423 
0424 if SideInfo == 1
0425   for z=1:r
0426     for w=1:r
0427       F_rfc_z{z,w} = c_m*F_rfc_z{z,w}; 
0428     end
0429   end
0430 end
0431 
0432 T(9,:)=clock;
0433 
0434

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