RFCDEMO1 Demo for switching AR(1)-processes. [F_RFC] = rfcdemo1; [F_RFC] = rfcdemo1(demoNr); [F_RFC] = rfcdemo1(demoNr,P); [F_RFC] = rfcdemo1(demoNr,P,A,m,s2,param); Input: demoNr = 1 / 2 / 3 / 0 Default: 1. 0 = Define your own process. P = Transition matrix for regime process. [rxr] A = AR coefficients. [rx2] m = Mean koefficients. [rx1] s2 = Innovation variances. [rx1] param = Parameter vector, [a b n], defines discretization. Output: F_RFC = Calculated rainflow intensity [nxn] Demo for switching AR(1)-processes, by using switching Markov chains. See Examples 3.1-3.3 in PhD-thesis: P. Johannesson (1999): Rainflow Analysis of Switching Markov Loads. [http://www.maths.lth.se/matstat/staff/pj/PhD/] Example: rfcdemo1; rfcdemo1(2); P=[0.98 0.02; 0.05 0.95]; rfcdemo1(1,P); See also rfcdemo2, fatigue
Computes the spectral density for an AR- or ARMA-model. | |
Calculates the number of upcrossings from a cycle count | |
Plots a cycle matrix, e.g. a rainflow matrix. | |
Plots cycles as points together with isolines of a cycle matrix. | |
Extracts turning points from data, | |
Discretizes an AR(1) process. | |
plots a Hidden Markov Model. | |
Calculates discrete levels given the parameter matrix. | |
Calculates the rainflow matrix/intensity for a Markov chain. | |
Calculates the stationary distribution for a Markov chain. | |
Simulates a switching ARMA-process. | |
Calculates the rainflow matrix/intensity for a switching Markov chain. | |
Finds the rainflow cycles from the sequence of turning points. | |
Root directory of WAFO installation. | |
Add directory to search path. | |
Control axis scaling and appearance. | |
Delete file or graphics object. | |
Flush pending graphics events. | |
Display message and abort function. | |
Check if variables or functions are defined. | |
Create figure window. | |
Directory separator for this platform. | |
Get handle to current figure. | |
Grid lines. | |
Hold current graph. | |
Convert number to string. (Fast version) | |
Wait for user response. | |
Linear plot. | |
Remove directory from search path. | |
Set object properties. | |
Square wave generation. | |
Stairstep plot. | |
Create axes in tiled positions. | |
Graph title. | |
Display warning message; disable or enable warning messages. | |
X-axis label. | |
Y-axis label. |
001 function [F_RFC] = refdemo1(demoNr,P,A,m,s2,param) 002 003 % RFCDEMO1 Demo for switching AR(1)-processes. 004 % 005 % [F_RFC] = rfcdemo1; 006 % [F_RFC] = rfcdemo1(demoNr); 007 % [F_RFC] = rfcdemo1(demoNr,P); 008 % [F_RFC] = rfcdemo1(demoNr,P,A,m,s2,param); 009 % 010 % Input: 011 % demoNr = 1 / 2 / 3 / 0 012 % Default: 1. 0 = Define your own process. 013 % P = Transition matrix for regime process. [rxr] 014 % A = AR coefficients. [rx2] 015 % m = Mean koefficients. [rx1] 016 % s2 = Innovation variances. [rx1] 017 % param = Parameter vector, [a b n], defines discretization. 018 % 019 % Output: 020 % F_RFC = Calculated rainflow intensity [nxn] 021 % 022 % Demo for switching AR(1)-processes, by using 023 % switching Markov chains. 024 % See Examples 3.1-3.3 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 % rfcdemo1; 031 % rfcdemo1(2); 032 % P=[0.98 0.02; 0.05 0.95]; rfcdemo1(1,P); 033 % 034 % See also rfcdemo2, fatigue 035 036 % Tested on Matlab 5.3 037 % 038 % History: 039 % Created by PJ (Pär Johannesson) 1997 040 % Updated by PJ 19-May-2000 041 % updated for WAFO 042 % Updated by PJ 07-Jul-2005 043 % warning off 044 045 warning_state=warning; % Get current warning-settings 046 warning off 047 048 fprintf(1,'Welcome to demo 1 for "Rainflow Cycles for Switching Processes"!\n'); 049 fprintf(1,'It demonstrates calculation of the rainflow intensity\n'); 050 fprintf(1,'for a switching AR(1)-process.\n'); 051 fprintf(1,'Type "help rfcdemo1" for further information\n\n'); 052 053 % Check input arguments 054 055 ni = nargin; 056 no = nargout; 057 error(nargchk(0,6,ni)); 058 059 Zstr='123456789'; 060 061 % Add path to rfcdemo1 062 demoDir = [waforoot filesep 'wdemos' filesep 'rfcdemo1']; 063 addpath(demoDir); 064 065 % Delete all figure windows 066 067 slut=0; 068 h=gcf; 069 while ~slut 070 h_prev=h; 071 delete(h); 072 h=gcf; 073 if h==h_prev, slut=1; end 074 end 075 076 077 if ni == 0, demoNr=1; end 078 079 % Define switching AR(1)-process 080 081 if demoNr == 1 082 p1=0.02; p2=0.01; 083 P0 = [1-p1 p1; p2 1-p2]; 084 A = [1 -0.5; 1 0.5]; 085 m = 2*[-0.5 1.5]'; 086 s2 = [1 1]'; 087 param = [-8 8 64]; 088 elseif demoNr == 2 089 p1=0.01; p2=0.02; 090 P0 = [1-p1 p1; p2 1-p2]; 091 A = [1 -0.5; 1 0.5]; 092 m = [0 0]'; 093 s2 = [4 1]'; 094 param = [-10 10 64]; 095 elseif demoNr == 3 096 p12=0.02; p13=0.005; 097 p21=0.01; p23=0.01; 098 p31=0.005; p32=0.02; 099 P0 = [1-p12-p13 p12 p13; p21 1-p21-p23 p23; p31 p32 1-p31-p32]; 100 A = [1 -0.5; 1 -0.3; 1 0.5]; 101 m = 2*[-0.5 0 1.5]'; 102 s2 = [1 1 1.44]'; 103 param = [-8 8 64]; 104 elseif demoNr == 0 105 if ni < 6 106 error('demoNr=0: Too few input arguments.'); 107 end 108 end 109 110 if ~exist('P') 111 P=P0; 112 end 113 114 % Initiation 115 116 u=levels(param); % Discretization levels 117 delta = u(2)-u(1); % Discretization step 118 r=length(P); % Number of regime states 119 n=param(3); % Number of discretization levels 120 C=ones(r,1); % 121 122 % Calculate statistics: mean , std, proportion, mean time 123 124 mX = m./sum(A')'; 125 sX = sqrt(s2./(1-A(:,2).^2)); 126 statP = mc2stat(P); 127 tau = 1./(1-diag(P))'; 128 129 % Parameters 130 131 fprintf('Parameters\n'); 132 fprintf(' z a(z) m(z) s(z)\n'); 133 fprintf('%2d %8g %8g %8g\n',[(1:r)' A(:,2) m sqrt(s2)]') 134 P 135 fprintf('\n') 136 137 % Statistics 138 fprintf('Statistics\n'); 139 fprintf('\n z m_X(z) s_X(z) ro(z) tau(z)\n'); 140 fprintf('%2d %8g %8g %8g %8g\n',[(1:r)' mX sX statP' tau']') 141 fprintf('\n') 142 143 % Calculate spectral densities 144 145 fprintf(1,' -- Calculate spectral densities.\n'); 146 147 set(gcf,'Name','Power spectra') 148 149 for i=1:r 150 Si = armaspec(1,A(i,:),s2(i,:),100); 151 subplot(r,1,i) 152 plot(Si(:,1),Si(:,2)), grid 153 ylabel(['S(f), regime ' Zstr(i)]) 154 end 155 xlabel('frequency, f') 156 157 drawnow, disp('Hit any key to continue.'); pause; 158 159 % Simulate MARX(2,1)-process 160 161 fprintf(1,' -- Simulating.\n'); 162 163 T=10000; % Length of simulation 164 Tinit=500; % Length of initiation (to remove transients) 165 166 [x,z,e] = sarmasim(C,A,m,s2,P,T,Tinit); 167 168 figure 169 set(gcf,'Name','Switching AR(1)-process') 170 171 I = 1:501; 172 hmmplot(x(I),z(I),I-1,[1 r],'','Switching Process, X(k)',1); 173 ylabel('regime, Z(k)'),xlabel('time, k') 174 175 drawnow, disp('Hit any key to continue.'); pause; 176 177 % 178 % Discretize AR(1)-processes 179 % Approximate with Markov chain 180 % 181 182 ud = u; ud(1)=u(1)-0.25*(u(n)-u(1)); ud(n)=u(n)+0.25*(u(n)-u(1)); 183 184 figure 185 set(gcf,'Name','Markov chain for Discretized AR(1)-processes') 186 187 if demoNr == 1 | demoNr ==2 | demoNr ==3 188 fprintf(1,' -- Loading discretized AR-processes.\n'); 189 190 eval(['load rfcdem1' Zstr(demoNr)]); % load uu,Q1,Q2,Q3 191 for i=1:r 192 subplot(1,r,i) 193 eval(['Q{i}=Q' num2str(i) ';']); 194 cmatplot(uu,uu,Q{i},3), axis square 195 title(['Regime ' Zstr(i)]) 196 end; 197 drawnow 198 else 199 200 fprintf(1,' -- Discretizing AR-processes. (This can take some time!)\n'); 201 for i=1:r 202 fprintf(1,' -- Discretizing process %d.\n',i); 203 a1i = A(i,2); 204 si = sqrt(s2(i)); 205 mi = mX(i); 206 [Qi,uu] = discar(ud,a1i,mi,si); 207 Q{i}=Qi; 208 %eval(['Q{i}' num2str(i) '=Qi;']); 209 subplot(r,1,i) 210 cmatplot(uu,uu,Qi,3), axis square 211 title(['Regime ' Zstr(i)]) 212 drawnow 213 end 214 end 215 216 drawnow, disp('Hit any key to continue.'); pause; 217 218 % 219 % Compute RFC counting distribution for SMC 220 % 221 222 fprintf(1,' -- Calculating rainflow intensity.\n'); 223 224 Qstr=''; 225 for i=1:r 226 Qstr = [Qstr ',Q' Zstr(i)]; 227 end 228 229 [F_RFC,mu_RFC] = smc2rfm(P,Q,[2 delta]); % definition 2 230 %eval( ['[F_RFC,mu_RFC] = smc2rfc(P,[2 delta]' Qstr ');'] ); % definition 2 231 232 for i=1:r 233 %eval( ['[F_RFC' Zstr(i) ',mu_RFC' Zstr(i) '] = mc2rfc(Q' Zstr(i) ',[2 delta]);'] ); % Regime i 234 [F_RFC0{i},mu_RFC0{i}] = mc2rfm(Q{i},[2 delta]); % Regime i 235 end 236 237 figure 238 239 cmatplot(u,u,F_RFC,2); 240 set(gcf,'Name','Rainflow intensity for switching AR(1)-process') 241 title('3D-plot') 242 243 drawnow, disp('Hit any key to continue.'); pause; 244 245 figure 246 247 % Calculate Cycles and crossings 248 249 fprintf(1,' -- Cycles and crossings for simulated load.\n'); 250 251 TP = dat2tp([(1:T)' x]); 252 RFC = tp2rfc(TP); 253 cross = cc2lc(RFC); 254 255 cocc(param,RFC,F_RFC) 256 set(gcf,'Name','Rainflow intensity for switching AR(1)-process') 257 title('Iso-lines compared with observed cycles') 258 xlabel('min'), ylabel('Max') 259 260 drawnow, disp('Hit any key to continue.'); pause; 261 262 % Plot Crossing intensity 263 264 figure 265 set(gcf,'Name','Crossing intensity') 266 267 subplot(1,2,1) 268 plot(u,diag(mu_RFC),'g--'), hold on 269 stairs(cross(:,1),cross(:,2)/T), hold off 270 title('Calculated vs. observed') 271 xlabel('level, u'); ylabel('Crossing intensity, mu(u)'); 272 273 subplot(1,2,2) 274 plot(u,diag(mu_RFC),'g'), hold on 275 cross_sum=zeros(size(diag(mu_RFC))); 276 for i=1:r 277 crossi=diag(mu_RFC0{i}); 278 plot(u,statP(i)*crossi,'r-.') 279 cross_sum=cross_sum+statP(i)*crossi; 280 end 281 plot(u,cross_sum,'y--'), hold off 282 title('Calculated vs. individual') 283 xlabel('level, u'); ylabel('Crossing intensity, mu(u)'); 284 285 fprintf(1,'\nThank you for running this demo!\n'); 286 287 % Remove path to rfcdemo1 288 demoDir = [waforoot filesep 'wdemos' filesep 'rfcdemo1']; 289 rmpath(demoDir); 290 291 warning(warning_state); 292 293
Comments or corrections to the WAFO group