Home > wafo > markov > smctp2rfm.m

smctp2rfm

PURPOSE ^

Calculates the rainflow matrix for a SMCTP.

SYNOPSIS ^

[F_rfc,mu_rfc,T] = smctp2rfm(P,F,c_m)

DESCRIPTION ^

  SMCTP2RFM  Calculates the rainflow matrix for a SMCTP.
 
  CALL: [Frfc,mu_rfc] = smctp2rfm(P,F,c_m);
 
  Frfc   = Rainflow Matrix                            [n,n]
  mu_rfc = Rainflow Counting intensity                [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)
  If a matrix F{i,2}=[], then the process will
  be assumed to be time-reversible.
 
  Calculates the rainflow matrix for a switching process
    with a Markov chain of turning points within each regime.
 
  Example: (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 = smctp2rfm(P,{F1 F1'; F2 F2'});
    cmatplot(u,u,Frfc)
 
  See also  mctp2rfm, smctp2arfm

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

001 function [F_rfc,mu_rfc,T] = smctp2rfm(P,F,c_m)
002 % SMCTP2RFM  Calculates the rainflow matrix for a SMCTP.
003 %
004 % CALL: [Frfc,mu_rfc] = smctp2rfm(P,F,c_m);
005 %
006 % Frfc   = Rainflow Matrix                            [n,n]
007 % mu_rfc = Rainflow Counting intensity                [n,n]
008 %
009 % P      = Transition matrix for regime process.      [r,r]
010 % F      = Cell array of min-Max and Max-min matrices {r,2}
011 % F{i,1} = min-Max matrix, process i                  [n,n]
012 % F{i,2} = Max-min matrix, process i                  [n,n]
013 % c_m    = Intensity of local minima, switching proc. [1,1]
014 %          (Default: 1)
015 % If a matrix F{i,2}=[], then the process will
016 % be assumed to be time-reversible.
017 %
018 % Calculates the rainflow matrix for a switching process
019 %   with a Markov chain of turning points within each regime.
020 %
021 % Example: (Example 4.1 in PhD thesis)
022 %   P = [0.9 0.1; 0.05 0.95];
023 %   param = [-1 1 32]; u = levels(param);
024 %   F1 = mktestmat(param,[-0.4 -0.3],0.15,1);
025 %   F2 = mktestmat(param,[0.3 0.4],0.15,1);
026 %   Frfc = smctp2rfm(P,{F1 F1'; F2 F2'});
027 %   cmatplot(u,u,Frfc)
028 %
029 % See also  mctp2rfm, smctp2arfm
030   
031 % References
032 %  
033 %  P. Johannesson (1999):
034 %  Rainflow Analysis of Switching Markov Loads.
035 %  PhD thesis, Mathematical Statistics, Centre for Mathematical Sciences,
036 %  Lund Institute of Technology.
037 %  
038 %  P. Johannesson (1998):
039 %  Rainflow Cycles for Switching Processes with Markov Structure.
040 %  Probability in the Engineering and Informational Sciences, 
041 %  Vol. 12, No. 2, pp. 143-175.
042   
043 % Tested  on Matlab  5.3
044 %
045 % History:
046 % Revised by PJ  23-Nov-1999
047 %   updated for WAFO
048 % Created by PJ (Pär Johannesson) 1997
049 %   Copyright (c) 1997-1998 by Pär Johannesson
050 %   Toolbox: Rainflow Cycles for Switching Processes V.1.1, 22-Jan-1998
051 
052 
053 % Check input arguments
054 
055 ni = nargin;
056 no = nargout;
057 error(nargchk(2,3,ni));
058 
059 if ni < 3
060   c_m=[];
061 end
062 
063 if isempty(c_m)
064   c_m=1;
065 end
066 
067 % Define 
068 
069 Zstr = '123456789'; Zstr=Zstr';
070 
071 T(1,:)=clock;
072 
073 r = length(P);   % Number of regime states
074 
075 n = length(F{1,1});  % Number of levels
076 
077 % Check that the rowsums of P are equal to 1
078 
079 sumP = sum(P');
080 if sum(sumP == 1) ~= r
081   warning(['Rowsums of P not equal to 1. Renormalizing!']);
082   for i = 1:length(P)
083     P(i,:) = P(i,:)/sumP(i);
084   end
085 end
086 
087 T(2,:)=clock;
088 
089 % Normalize the rowsums of F{1,1},...,F{r,1} to 1
090 %  ==>  Q{1,1},...,Q{r,1}
091 
092 for i = 1:r
093   QQ{i,1} = triu(F{i,1},1); % Set zeros below diagonal and diagonal
094   % Set negative elements to zero
095   [I,J] = find(QQ{i,1}<0);
096   if length(I) ~= 0
097     warning(['Negative elements in Q' Zstr(i) '. Setting to zero!']);
098     for k = 1:length(I)
099       QQ{i,1}(I(k),J(k)) = 0;
100     end
101   end
102 
103   sumQQi = sum(QQ{i,1}');
104   % Normalize rowsums
105   if sum(sumQQi == 1) ~= length(QQ{i,1})
106     %disp(['Warning: Rowsums of Q' Zstr(i) ' not equal to 1. Renormalizing!']);
107     for j = 1:n-1
108       if sumQQi(j)~=0, QQ{i,1}(j,:) = QQ{i,1}(j,:)/sumQQi(j); end
109     end
110   end
111 end
112 
113 T(3,:)=clock;
114 
115 % Normalize the rowsums of F{1,2},...,F{r,2} to 1
116 %  ==>  Q{1,2},...,Q{r,2}
117 
118 % Normalize the rowsums of Fh1,...,Fhr to 1  ==>  Q1,...,Qr
119 
120 for i = 1:r
121   if isempty(F{i,2})        % Time-reversible
122     QQ{i,2} = F{i,1}';
123   else                   % Fhi is given
124     QQ{i,2} = F{i,2}; 
125   end
126     
127   QQ{i,2} = tril(QQ{i,2},-1); % Set zeros above diagonal and diagonal
128   % Set negative elements to zero
129   [I,J] = find(QQ{i,2}<0);
130   if length(I) ~= 0
131     warning(['Negative elements in Qh' Zstr(i) '. Setting to zero!']);
132     for k = 1:length(I)
133       QQ{i,2}(I(k),J(k)) = 0;
134     end
135   end
136 
137   sumQQi = sum(QQ{i,2}');
138   if sum(sumQQi == 1) ~= length(QQ{i,2})
139     %disp(['Warning: Rowsums of Qh' Zstr(i) ' not equal to 1. Renormalizing!']);
140     for j = 2:n
141       if sumQQi(j)~=0, QQ{i,2}(j,:) = QQ{i,2}(j,:)/sumQQi(j); end
142     end
143   end
144 end
145 
146 T(4,:)=clock;
147 
148 % Make the transition matrix Q for the joint min-Max process
149 
150 Q = zeros(n*r,n*r);
151 I = 0:r:(n-1)*r;
152 for z=1:r
153   Q0 = kron(QQ{z,1},P);
154   Q(I+z,:) = Q0(I+z,:);
155 end
156 
157 T(5,:)=clock;
158 
159 % Make the transition matrix Qh for the joint Max-min process
160 
161 Qh = zeros(n*r,n*r);
162 I = 0:r:(n-1)*r;
163 for z=1:r
164   Q0 = kron(QQ{z,2},P);
165   Qh(I+z,:) = Q0(I+z,:);
166 end
167 
168 T(6,:)=clock;
169 
170 % Stationary distribution (=ro) of local minima with transition matrix
171 % Qt = Q*Qh = "Transition matrix for min-to-min"
172 
173 Qt = Q*Qh;
174 ro = mc2stat(Qt(1:r*(n-1),1:r*(n-1)));  % Stationary distr., row vector  
175 ro = [ro zeros(1,r)];  % Minimum can't reach the highest level
176 
177 T(7,:)=clock;
178 
179 % Calculate rainflow matrix 
180 
181 mu_rfc = zeros(n,n);
182 EYE = eye(n*r,n*r);
183 Qcumsum = fliplr(cumsum(fliplr(Q)')');
184 
185 %fprintf(1,'Calculating row ');
186 for i=2:n
187   %fprintf(1,'-%1d',i);
188   for j=i-1:min([i n-1])
189     q = Qcumsum(:,r*j+1);
190 
191     d = q(1:r*(i-1));                           % 1:i-1
192     Ro = ro(1:r*(i-1));                         % 1:i-1
193 
194     mu_rfc(i,j) = Ro*d;
195   end
196   for j=i+1:n-1
197     q = Qcumsum(:,r*j+1);
198 
199     A = Q(r*(i-1)+1:r*(j-1),r*(i+1-1)+1:r*j);   % i:j-1, i+1:j
200     Ah = Qh(r*(i+1-1)+1:r*j,r*(i-1)+1:r*(j-1)); % i+1:j, i:j-1
201     C = Q(1:r*(i-1),r*(i+1-1)+1:r*j);           % 1:i-1, i+1:j
202     Eye = EYE(1:r*(j-i),1:r*(j-i));             % eye (size of A*Ah)
203     d = q(1:r*(i-1));                           % 1:i-1
204     e = q(r*(i-1)+1:r*(j-1));                   % i:j-1
205     Ro = ro(1:r*(i-1));                         % 1:i-1
206 
207     mu_rfc(i,j) = Ro*(d+C*Ah*inv(Eye-A*Ah)*e);
208   end
209 end
210 %fprintf(1,'\n');
211 
212 
213 % Calculation of the intensity of local minima for the swithching process
214 
215 mu_rfc = c_m*mu_rfc;
216 
217 T(8,:)=clock;
218 
219 % Convert: mu_rfc  -->  F_rfc
220 
221 F_rfc = zeros(n,n);
222 for i= 1:n-1
223   for j= i+1:n
224     F_rfc(i,j) = mu_rfc(i+1,j-1) - mu_rfc(i,j-1) - mu_rfc(i+1,j) + mu_rfc(i,j);
225   end
226 end
227 
228 I = F_rfc<0;
229 if sum(sum(I)) ~= 0
230   warning(['Negative elements in calculated rainflow matrix F_rfc. Setting to zero!']);
231   F_rfc(I) = 0;
232 end
233 
234 
235 if no>1
236   mu_rfc = cmat2nt(F_rfc); % Calculate rainflow counting intensity
237 end
238 
239 
240 T(9,:)=clock;
241 
242

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