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)

## CROSS-REFERENCE INFORMATION

This function calls:
 cmat2nt Calculates a counting distribution from a cycle 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. kron Kronecker tensor product. tril Extract lower triangular part. triu Extract upper triangular part. warning Display warning message; disable or enable warning messages.
This function is called by:
 f_smctp Auxiliary function used by ESTSMCTP itmkurs_lab2 Script to computer exercises 2 itmkurs_lab4 Script to computer exercises 4 mctp2rfm Calculates the rainflow matrix for a MCTP. rfcdemo2 Rainflow matrix for Switching Markov Chains of Turning Points. test_markov Quick test of the routines in module 'markov'

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