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 
      Switching Markov Loads. 
      [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); 
  
  See also  rfcdemo1, fatigue

CROSS-REFERENCE INFORMATION ^

This function calls: 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 
026 %     Switching Markov Loads. 
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 % 
038 % See also  rfcdemo1, fatigue 
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  
158 % Simulate a load 
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