Home > wafo > wdemos > rfcdemo2.m

# rfcdemo2

## PURPOSE

Rainflow matrix for Switching Markov Chains of Turning Points.

## SYNOPSIS

[F_RFC] = refdemo2(demoNr,P,param,x0,s,lam)

## DESCRIPTION

```  RFCDEMO2 Rainflow matrix for Switching Markov Chains of Turning Points.

[F_RFC] = rfcdemo2;
[F_RFC] = rfcdemo2(P);
[F_RFC] = rfcdemo2(demoNr,P,param,x0,s,lam);

Input:
demoNr = 1:  Example 4.1  (Default: 1)
2:  Example 4.2
11: Another example
P      = Transition matrix for regime process.  [rxr]
param  = Parameter vector, [a b n], defines discretization.
x0     = Center of ellipse [min Max]            [rx2]
s      = Standard deviation (0<s<infinity)      [rx1]
lam    = Orientation of ellipse (0<lam<infinity)[rx1]

Output:
F_RFC  = Calculated rainflow intensity          [nxn]

The demo presents calculation of the rainflow matrix for
a Switching Markov Chains of Turning Points (SMCTP).
See Examples 4.1 and 4.2 in PhD-thesis:
P. Johannesson (1999): Rainflow Analysis of
[http://www.maths.lth.se/matstat/staff/pj/PhD/]

Example:
rfcdemo2;
rfcdemo2(1);
P=[0.98 0.02; 0.05 0.95];
rfcdemo2(1,P);
P=[0.98 0.02; 0.05 0.95]; x0=[-0.4 -0.3; 0.3 0.4]; param=[-1 1 32];
s=[0.1 0.15]'; lam=[0.7 1.5]';
rfcdemo2(1,P,param,x0,s,lam);

## CROSS-REFERENCE INFORMATION

This function calls:
 cc2cmat Calculates the cycle count matrix from a cycle count. cc2lc Calculates the number of upcrossings from a cycle count ccplot Plots a cycle count as a point process in the plane. cmat2dmat Computes the (Palmgren-Miner) damage matrix from a cycle matrix. cmatplot Plots a cycle matrix, e.g. a rainflow matrix. dat2tp Extracts turning points from data, hmmplot plots a Hidden Markov Model. levels Calculates discrete levels given the parameter matrix. mc2stat Calculates the stationary distribution for a Markov chain. mktestmat Makes test matrices for min-max (and max-min) matrices. smctp2rfm Calculates the rainflow matrix for a SMCTP. smctpsim Simulates a switching Markov chain of turning points, smoothcmat Smooth a cycle matrix using (adaptive) kernel smoothing tp2rfc Finds the rainflow cycles from the sequence of turning points. axis Control axis scaling and appearance. cell Create cell array. delete Delete file or graphics object. drawnow Flush pending graphics events. error Display message and abort function. figure Create figure window. gcf Get handle to current figure. grid Grid lines. hold Hold current graph. pause Wait for user response. plot Linear plot. set Set object properties. subplot Create axes in tiled positions. title Graph title. warning Display warning message; disable or enable warning messages. xlabel X-axis label. ylabel Y-axis label.
This function is called by:

## SOURCE CODE

```001 function [F_RFC] = refdemo2(demoNr,P,param,x0,s,lam)
002
003 % RFCDEMO2 Rainflow matrix for Switching Markov Chains of Turning Points.
004 %
005 % [F_RFC] = rfcdemo2;
006 % [F_RFC] = rfcdemo2(P);
007 % [F_RFC] = rfcdemo2(demoNr,P,param,x0,s,lam);
008 %
009 % Input:
010 % demoNr = 1:  Example 4.1  (Default: 1)
011 %          2:  Example 4.2
012 %          11: Another example
013 % P      = Transition matrix for regime process.  [rxr]
014 % param  = Parameter vector, [a b n], defines discretization.
015 % x0     = Center of ellipse [min Max]            [rx2]
016 % s      = Standard deviation (0<s<infinity)      [rx1]
017 % lam    = Orientation of ellipse (0<lam<infinity)[rx1]
018 %
019 % Output:
020 % F_RFC  = Calculated rainflow intensity          [nxn]
021 %
022 %   The demo presents calculation of the rainflow matrix for
023 %   a Switching Markov Chains of Turning Points (SMCTP).
024 %   See Examples 4.1 and 4.2 in PhD-thesis:
025 %     P. Johannesson (1999): Rainflow Analysis of
027 %     [http://www.maths.lth.se/matstat/staff/pj/PhD/]
028 %
029 % Example:
030 %   rfcdemo2;
031 %   rfcdemo2(1);
032 %   P=[0.98 0.02; 0.05 0.95];
033 %   rfcdemo2(1,P);
034 %   P=[0.98 0.02; 0.05 0.95]; x0=[-0.4 -0.3; 0.3 0.4]; param=[-1 1 32];
035 %   s=[0.1 0.15]'; lam=[0.7 1.5]';
036 %   rfcdemo2(1,P,param,x0,s,lam);
037 %
039
040 % Tested  on Matlab  5.3
041 %
042 % History:
043 % Updated by PJ 18-May-2000
044 %   updated for WAFO
045 % Created by PJ (Pär Johannesson) 1997
046 % Updated by PJ 07-Jul-2005
047 %   warning off
048
049 warning_state=warning;  % Get current warning-settings
050 warning off
051
052 % Check input arguments
053
054 ni = nargin;
055 no = nargout;
056 error(nargchk(0,6,ni));
057
058
059 fprintf(1,'Welcome to demo 2 for "Rainflow Cycles for Switching Processes"!\n');
060 fprintf(1,'It demonstrates calculation of the rainflow matrix\n');
061 fprintf(1,'for a Switching Markov Chain of Turning Points.\n');
062 fprintf(1,'Type "help rfcdemo2" for further information\n\n');
063
064 Zstr='123456789';
065
066 % Delete all figure windows
067
068 slut=0;
069 h=gcf;
070 while ~slut
071   h_prev=h;
072   delete(h);
073   h=gcf;
074   if h==h_prev, slut=1; end
075 end
076
077 % Check input-arguments
078
079 if ni<1,        demoNr=1; end
080 if isempty(demoNr), demoNr=1; end
081 if ni<2,        P=[];     end
082 if ni<3,        param=[]; end
083 if ni<4,        x0=[];    end
084 if ni<5,        s=[];     end
085 if ni<6,        lam=[];   end
086
087 if demoNr == 1
088   p1=0.10; p2=0.05;
089   P0=[1-p1 p1; p2 1-p2];
090   param0 = [-1 1 32];
091   x00 = [-0.4 -0.3; 0.3 0.4];
092
093   s0 = [0.15 0.15]';
094   lam0 = [1.0 1.0]';
095 elseif demoNr == 2
096   p1=0.10; p2=0.05;
097   P0=[1-p1 p1; p2 1-p2];
098   param0 = [-1 1 32];
099   x00 = [-0.1 0.1; 0.0 0.0];
100
101   s0 = [0.28 0.12]';
102   lam0 = [0.5 2]';
103 elseif demoNr == 11
104   p1=0.04; p2=0.01;
105   P0=[1-p1 p1; p2 1-p2];
106   param0 = [-1 1 32];
107   x00 = [0.2 0.3; -0.1 -0.1];
108   s0 = [0.18 0.1]';
109   lam0 = [.5 2]';
110 end
111
112 if isempty(P),     P=P0;         end
113 if isempty(param), param=param0; end
114 if isempty(x0),    x0=x00;       end
115 if isempty(s),     s=s0;         end
116 if isempty(lam),   lam=lam0;     end
117
118 % Define min-Max and Max-min matrices
119
120 % Initiation
121
122 u=levels(param);    % Discretization levels
123 r=length(P);        % Number of regime states
124 n=param(3);         % Number of discretization levels
125
126 paramD = [1 n n];
127
128 F=cell(r,2);
129 for i=1:r
130   [F{i,1},F{i,2}] = mktestmat(param,x0(i,:),s(i),lam(i));
131 end
132
133
134 % Calculate statistics: mean , std, proportion, mean time
135
136 %mX = m./sum(A')';
137 %sX = sqrt(s2./(1-A(:,2).^2));
138 statP = mc2stat(P);
139 tau = 1./(1-diag(P))';
140
141
142 % Parameters
143
144 fprintf('Parameters\n');
145 %fprintf(' z      a(z)      m(z)      s(z)\n');
146 %fprintf('%2d  %8g  %8g  %8g\n',[(1:r)' A(:,2) m sqrt(s2)]')
147 P
148
149 % Statistics
150 fprintf('Statistics\n');
151 fprintf('\n z    ro(z)    tau(z)\n');
152 fprintf('%2d  %8g  %8g\n',[(1:r)' statP' tau']')
153 fprintf('\n')
154
155 %fprintf('\n z    m_X(z)    s_X(z)     ro(z)    tau(z)\n');
156 %fprintf('%2d  %8g  %8g  %8g  %8g\n',[(1:r)' mX sX statP' tau']')
157
159
160 fprintf(1,' -- Simulating.\n');
161
162 T = 5000;           % Number of TP
163
164 [xD,z] = smctpsim(P,F,T);
165 x = u(xD)';
166
167 I = 1:501;
168 hmmplot(x(I),z(I),I-1,[1 r],'','Switching Process, X(k)',1);
169 ylabel('regime, Z(k)'),xlabel('time, k')
170
171 set(gcf,'Name','Switching MCTP')
172
173 drawnow, disp('Hit any key to continue.'); pause;
174
175 % Calculate the rainflow matrix
176
177 fprintf(1,' -- Calculating rainflow matrix.\n');
178
179 F_RFCsum = zeros(n,n);
180 for i=1:r
181   [F_RFC1{i},mu_RFC1{i}] = smctp2rfm(1,F(i,1:2));
182   F_RFCsum = F_RFCsum + statP(i)*F_RFC1{i};
183 end
184
185 figure
186 set(gcf,'Name','Calculated Rainflow matrices')
187 v = [param(1:2) param(1:2) 0 max(max(F_RFCsum))];
188 for i=1:r
189   subplot(1,r,i),
190   cmatplot(u,u,statP(i)*F_RFC1{i},1);
191   grid,axis('square'),axis(v)
192   title(['Regime ' Zstr(i)])
193 end
194
195 [F_RFC,mu_RFC] = smctp2rfm(P,F);
196
197 figure
198 set(gcf,'Name','Calculated Rainflow matrix')
199
200 v = [param(1:2) param(1:2) 0 max(max(F_RFCsum))];
201 subplot(1,2,1)
202 cmatplot(u,u,F_RFC,1),grid
203 axis('square'),axis(v)
204 title('Mixed RFM')
205 subplot(1,2,2)
206 cmatplot(u,u,F_RFCsum,1),grid
207 axis('square'),axis(v)
208 title('Sum of RFMs')
209
210 drawnow, disp('Hit any key to continue.'); pause;
211
212 figure
213
214 % Calculate Cycles and crossings
215
216 fprintf(1,' -- Cycles and crossings for simulated load.\n');
217
218 TP = dat2tp([(1:T)' x]);
219 RFC = tp2rfc(TP);
220 ccplot(RFC)
221 cross = cc2lc(RFC,1);   % Only upcrossings
222
223 % Observed Rainflow Matrix
224
225 TPD = dat2tp([(1:T)' xD]);
226 RFCD = tp2rfc(TPD);
227 fr = cc2cmat(paramD,RFCD);
228 F_RFCobs = fr;
229 %[F_RFCobs_smooth,h] = smthcmat(F_RFCobs);
230 [F_RFCobs_smooth,h] = smoothcmat(F_RFCobs,1,0.7);
231
232 set(gcf,'Name','Observed Rainflow Matrix')
233 subplot(1,2,1), cmatplot(u,u,F_RFCobs/T*2,1),grid
234 axis('square'),axis(v)
235 title('Observed RFM')
236 subplot(1,2,2), cmatplot(u,u,F_RFCobs_smooth/T*2,1), grid
237 axis('square'),axis(v)
238 title('Smoothed observed RFM')
239
240 drawnow, disp('Hit any key to continue.'); pause;
241
242
243 % Check crossing intensity
244
245 % Check Damage matrix
246
247 figure
248 set(gcf,'Name','Damage')
249 beta0 = 3;
250 Dmat = cmat2dmat(param,F_RFC,beta0);
251 Dmatsum = cmat2dmat(param,F_RFCsum,beta0);
252 subplot(1,2,1),cmatplot(u,u,Dmat),grid
253 axis('square')
254
255 v=axis;
256 title('Mixed RFM')
257 subplot(1,2,2),cmatplot(u,u,Dmatsum),grid
258 axis('square')
259
260 axis(v);
261 title('Sum of RFMs')
262
263
264 drawnow, disp('Hit any key to continue.'); pause;
265
266
267 % Compare crossing spectrum with simulated
268
269 figure
270 set(gcf,'Name','Crossing Intensity')
271
272 subplot(1,2,1)
273 plot(u,diag(mu_RFC),'g--'), hold on
274 %[XX,YY]=stairs(u,diag(mu_RFC));
275 %plot(XX,YY,'g'),hold on
276 %stairs(cross(:,1),cross(:,2)/T*2), hold off
277 plot(cross(:,1),cross(:,2)/T*2), hold off
278 xlabel('level, u'); ylabel('mu(u)');
279 v = axis; axis([param(1:2) 0 v(4)]);
280
281 subplot(1,2,2)
282 plot(u,diag(mu_RFC),'g'), hold on
283 plot(u,diag(mu_RFC1{1})*statP(1),'r-.')
284 plot(u,diag(mu_RFC1{2})*statP(2),'r-.')
285 plot(u,diag(mu_RFC1{1})*statP(1)+diag(mu_RFC1{2})*statP(2),'y--'), hold off
286 xlabel('level, u'); ylabel('mu(u)');
287 v = axis; axis([param(1:2) 0 v(4)]);
288
289 fprintf(1,'\nThank you for running this demo!\n');
290
291 warning(warning_state);
292
293
294```

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