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

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. deal Deal inputs to outputs. error Display message and abort function. inv Matrix inverse. kron Kronecker tensor product. num2str Convert number to string. (Fast version)
This function is called by:
 f_smctp Auxiliary function used by ESTSMCTP

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