Home > wafo > markov > smc2rfm.m

smc2rfm

PURPOSE ^

Calculates the rainflow matrix/intensity for a switching Markov chain.

SYNOPSIS ^

[F_rfc,mu_rfc] = smc2rfc(P,Qc,def)

DESCRIPTION ^

 SMC2RFM Calculates the rainflow matrix/intensity for a switching Markov chain. 
  
  [F_rfc,mu_rfc] = smc2rfc(P,Q,1); 
  [F_rfc,mu_rfc] = smc2rfc(P,Q,[2 h]); 
  
  F_rfc  = Rainflow matrix / Rainflow intensity     [NxN] 
  mu_rfc = Rainflow counting intensity              [NxN] 
  
  P      = Transition matrix for regime process     [rxr] 
  Q      = Cell array of transition matrices        {r,1} 
  Q{i}   = Transition matrix for Markov chain i     [nxn] 
  def    = Definition 1: Markov chain (default),     N=n 
                      2: Discretized Markov chain,   N=n+1 
  h      = Discretization step (ONLY Def 2!) 
  
  Calculates  
    (1) the rainflow matrix for a switching Markov chain OR 
    (2) the rainflow intensity for a discretized switching Markov chain. 
  
  Example:  
    param = [-1 1 32]; u = levels(param); 
    [F1,F2] = mktestmat(param,[-0.3 0.3],0.15,1,-Inf); 
    Frfc = smc2rfm(P,{F1;F2}); 
    cmatplot(u,u,Frfc) 
  
  See also  mc2rfm, smctp2rfm, smc2stat, cmatplot

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

001 function [F_rfc,mu_rfc] = smc2rfc(P,Qc,def) 
002 %SMC2RFM Calculates the rainflow matrix/intensity for a switching Markov chain. 
003 % 
004 % [F_rfc,mu_rfc] = smc2rfc(P,Q,1); 
005 % [F_rfc,mu_rfc] = smc2rfc(P,Q,[2 h]); 
006 % 
007 % F_rfc  = Rainflow matrix / Rainflow intensity     [NxN] 
008 % mu_rfc = Rainflow counting intensity              [NxN] 
009 % 
010 % P      = Transition matrix for regime process     [rxr] 
011 % Q      = Cell array of transition matrices        {r,1} 
012 % Q{i}   = Transition matrix for Markov chain i     [nxn] 
013 % def    = Definition 1: Markov chain (default),     N=n 
014 %                     2: Discretized Markov chain,   N=n+1 
015 % h      = Discretization step (ONLY Def 2!) 
016 % 
017 % Calculates  
018 %   (1) the rainflow matrix for a switching Markov chain OR 
019 %   (2) the rainflow intensity for a discretized switching Markov chain. 
020 % 
021 % Example:  
022 %   param = [-1 1 32]; u = levels(param); 
023 %   [F1,F2] = mktestmat(param,[-0.3 0.3],0.15,1,-Inf); 
024 %   Frfc = smc2rfm(P,{F1;F2}); 
025 %   cmatplot(u,u,Frfc) 
026 % 
027 % See also  mc2rfm, smctp2rfm, smc2stat, cmatplot 
028  
029 % References 
030 %   
031 %  P. Johannesson (1999): 
032 %  Rainflow Analysis of Switching Markov Loads. 
033 %  PhD thesis, Mathematical Statistics, Centre for Mathematical Sciences, 
034 %  Lund Institute of Technology. 
035 %   
036 %  P. Johannesson (1998): 
037 %  Rainflow Cycles for Switching Processes with Markov Structure. 
038 %  Probability in the Engineering and Informational Sciences,  
039 %  Vol. 12, No. 2, pp. 143-175. 
040  
041 % Tested  on Matlab  5.3 
042 % 
043 % History: 
044 % Revised by PJ 19-May-2000 
045 %   Changer disp(...) to warning(...) . 
046 % Revised by PJ  23-Nov-1999 
047 %   updated for WAFO 
048 % Created by PJ (Pär Johannesson) 1997 
049 %   Copyright (c) 1997 by Pär Johannesson 
050 %   Toolbox: Rainflow Cycles for Switching Processes V.1.0, 2-Oct-1997 
051  
052  
053 % Check input arguments 
054  
055 ni = nargin; 
056 no = nargout; 
057 error(nargchk(2,3,ni)); 
058  
059 if ni<3, def = []; end 
060 if isempty(def), def = 1; end 
061  
062 % Define  
063  
064 Zstr = '123456789'; 
065  
066 r = length(P);   % Number of regime states 
067 n = length(Qc{1});  % Number of levels 
068  
069 % Check that the rowsums of P are equal to 1 
070  
071 sumP = sum(P'); 
072 if sum(sumP == 1) ~= length(P) 
073   warning(['Rowsums of P not equal to 1. Renormalizing.']); 
074   for i = 1:length(P) 
075     P(i,:) = P(i,:)/sumP(i); 
076   end 
077 end 
078  
079 % Check that the rowsums of Qc{1},...,Qc{r} are equal to 1 
080  
081 for i = 1:r 
082   sumQi = sum(Qc{i}'); 
083   if sum(sumQi == 1) ~= length(Qc{i}) 
084     warning(['Rowsums of Q{' Zstr(i) '} not equal to 1. Renormalizing.']); 
085     for j = 1:length(Qc{i}) 
086       Qc{i}(j,:) = Qc{i}(j,:)/sumQi(j); 
087     end 
088   end 
089 end 
090  
091  
092 % Make the transition matrix Q for the joint MC (X_k,Z_k) 
093  
094 Q = zeros(n*r,n*r); 
095 I = 0:r:(n-1)*r; 
096 for z=1:r 
097   QQ = kron(Qc{z},P); 
098   Q(I+z,:) = QQ(I+z,:); 
099 end 
100  
101  
102 % Stationary distribution = ro of Q 
103  
104 ro = mc2stat(Q); 
105  
106 % Calculate rainflow matrix / rainflow intensity 
107  
108 if def(1) == 1  % Markov Chain 
109  
110   N = n; 
111   mu_rfc = zeros(N,N); 
112   EYE = eye(n*r,n*r); 
113   Qcumsum = fliplr(cumsum(fliplr(Q)')'); 
114  
115   for i=2:n 
116     for j=i-1:n-1 
117       q = Qcumsum(:,r*j+1); 
118  
119       A = Q(r*(i-1)+1:r*j,r*(i-1)+1:r*j); % i:j,   i:j 
120       C = Q(1:r*(i-1),r*(i-1)+1:r*j);     % 1:i-1, i:j 
121       Eye = EYE(1:r*(j-i+1),1:r*(j-i+1)); % eye (size of A) 
122       d = q(1:r*(i-1));                   % 1:i-1 
123       e = q(r*(i-1)+1:r*j);               % i:j 
124       Ro = ro(1:r*(i-1));                 % 1:i-1 
125  
126       if j == i-1  % isempty(A) 
127         mu_rfc(i,j) = Ro*d; 
128       else 
129         mu_rfc(i,j) = Ro*(d+C*inv(Eye-A)*e); 
130       end 
131     end 
132   end 
133  
134 elseif def(1) == 2  % Discretized Markov Process 
135  
136   N = n+1; 
137   mu_rfc = zeros(N,N); 
138   EYE = eye(n*r,n*r); 
139   Qcumsum = fliplr(cumsum(fliplr(Q)')'); 
140  
141   for i=2:N-1 
142     for j=i:N-1 
143       q = Qcumsum(:,r*(j-1)+1); 
144  
145       A = Q(r*(i-1)+1:r*(j-1),r*(i-1)+1:r*(j-1)); % i:j-1, i:j-1 
146       C = Q(1:r*(i-1),r*(i-1)+1:r*(j-1));         % 1:i-1, i:j-1 
147       Eye = EYE(1:r*(j-i),1:r*(j-i));             % eye (size of A) 
148       d = q(1:r*(i-1));                           % 1:i-1 
149       e = q(r*(i-1)+1:r*(j-1));                   % i-1:j-1 
150       Ro = ro(1:r*(i-1));                         % 1:i-1 
151  
152       if i == j  % isempty(A) 
153         mu_rfc(i,j) = Ro*d; 
154       else 
155         mu_rfc(i,j) = Ro*(d+C*inv(Eye-A)*e); 
156       end 
157     end 
158   end 
159  
160 end 
161  
162 % Fill the missing triangel in the mu_rfc matrix 
163 if def(1) == 1 
164  
165   F_rfc = zeros(n); 
166   for i = 1:n-1 
167     for j= i+1:n 
168       F_rfc(i,j) = mu_rfc(i+1,j-1)-mu_rfc(i,j-1)-mu_rfc(i+1,j)+mu_rfc(i,j); 
169     end 
170   end 
171    
172   mu_rfc = cmat2nt(F_rfc); 
173    
174 % Old version 
175 %  mu_rfc=flipud(mu_rfc'); % Convert to WAT matrix-format 
176 %    for k = 2:N-1 
177 %      for m = 1:N-k 
178 %        i = k+m; 
179 %        j = N+1+k-i; 
180 %        mu_rfc(i,j) = mu_rfc(i,j-1)+mu_rfc(i-1,j)-mu_rfc(i-1,j-1); 
181 %      end 
182 %    end 
183 %     
184 %  mu_rfc = flipud(mu_rfc)'; % Convert to PJ-def 
185 %   
186 %  F_rfc = nt2cmat(mu_rfc);    % Calculate Rainflow matrix 
187 %   
188 %%  F_rfc = fliplr(triu(fliplr(F_rfc),1)); % Zeros blow SW-NE diagonal 
189  
190 end 
191  
192 if def(1) == 2 
193  
194   mu_rfc=flipud(mu_rfc'); % Convert to WAT matrix-format 
195   for k = 1:N-1 
196     for m = 1:N-k 
197       i = k+m; 
198       j = N+1+k-i; 
199       mu_rfc(i,j) = mu_rfc(i,j-1)+mu_rfc(i-1,j)-mu_rfc(i-1,j-1); 
200     end 
201   end 
202  
203   % Calculate rainflow intensity by numerical differentiation 
204  
205   h = def(2); 
206   F_rfc = zeros(N,N); 
207  
208   for i=2:N-1 
209     for j=2:N-i 
210       F_rfc(i,j) = (mu_rfc(i-1,j-1)-mu_rfc(i-1,j+1) -(mu_rfc(i+1,j-1)-mu_rfc(i+1,j+1))) / (2*h)^2; 
211     end 
212   end 
213    
214   F_rfc = flipud(F_rfc)';   % Convert to PJ-def 
215   mu_rfc = flipud(mu_rfc)'; % Convert to PJ-def 
216  
217 end 
218  
219  
220  
221

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