Home > wafo > wdemos > itmkurs > itmkurs_lab2.m

itmkurs_lab2

PURPOSE ^

Script to computer exercises 2

SYNOPSIS ^

This is a script file.

DESCRIPTION ^

 ITMKURS_LAB2 Script to computer exercises 2 
  
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
  lab2.m 
  Switching Markov Loads and Rainflow Analysis
 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

001 %ITMKURS_LAB2 Script to computer exercises 2 
002 % 
003 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
004 % lab2.m 
005 % Switching Markov Loads and Rainflow Analysis
006 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
007 
008 % 
009 % Section: Markov Chains of Turning Points
010 %
011 
012 % Model Definition
013 
014 n = 32; param = [-1 1 n];                % Define discretization
015 u = levels(param);                       % Discrete levels
016 G = mktestmat(param,[-0.2 0.2],0.15,1);  % min-max matrix
017 
018 % Simulation
019 
020 T = 5000;                % Length of simulation (number of TP)
021 xD = mctpsim({G []},T);  % Simulate
022 x = u(xD)';              % Change scale to levels -1,..,1
023 
024 t=1:200; plot(t,x(t))  
025 
026 % Calculating the Rainflow Matrix
027 
028 Grfc=mctp2rfm({G,[]});
029 
030 subplot(1,2,1),cmatplot(u,u,G),axis('square')
031 subplot(1,2,2),cmatplot(u,u,Grfc),axis('square')
032 
033 cmatplot(u,u,{G Grfc},3)
034 
035 FrfcObs = dtp2rfm(xD,n);
036 cmatplot(u,u,{FrfcObs Grfc*T/2})
037 
038 % Switching Markov Chains of Turning Points
039 
040 % Model Definition
041 
042 n=32; param = [-1 1 n];                   % Define discretization
043 u=levels(param);                          % Discrete levels
044 G1 = mktestmat(param,[-0.4 -0.3],0.15,1); % regime 1
045 G2 = mktestmat(param,[0.3 0.4],0.15,1);   % regime 2
046 
047 p1=0.10; p2=0.05;
048 P=[1-p1 p1; p2 1-p2]  % Transition matrix
049 statP=mc2stat(P)      % Stationary distribution
050 
051 % Simulation
052 
053 T=5000;   % Length of simulation (number of TP)
054 [xD,z] = smctpsim(P,{G1 []; G2 []},T); % Simulate
055 x=u(xD)'; % Change scale to levels -1,..,1
056 
057 t=1:400;
058 hmmplot(x(t),z(t),t,[1 2])         % Same colour
059 hmmplot(x(t),z(t),t,[1 2],'','',1) % Different colours
060 
061 % Calculating the Rainflow Matrix
062 
063 Grfc=smctp2rfm(P,{G1,[];G2,[]});
064 Grfc1=mctp2rfm({G1,[]});
065 Grfc2=mctp2rfm({G2,[]});
066 GrfcSum=statP(1)*Grfc1+statP(2)*Grfc2;
067 
068 cmatplot(u,u,{Grfc1 Grfc2; GrfcSum Grfc})
069 
070 % Model B (Optional)
071 
072 G1B = mktestmat(param,[-0.1 -0.1],0.28,0.5);  % regime 1
073 G2B = mktestmat(param,[0 0],0.12,2);          % regime 2
074 
075 [xDB,zB] = smctpsim(P,{G1B []; G2B []},T); % Simulate
076 xB=u(xDB)'; % Change scale to levels -1,..,1
077 hmmplot(xB(t),zB(t),t,[1 2],'','',1) % Different colours
078 
079 GrfcB=smctp2rfm(P,{G1B,[];G2B,[]});
080 Grfc1B=mctp2rfm({G1B,[]});
081 Grfc2B=mctp2rfm({G2B,[]});
082 GrfcSumB=statP(1)*Grfc1B+statP(2)*Grfc2B;
083 
084 % Observed Rainflow Matrix and Smoothing
085 
086 TP = dat2tp([(1:T)' xD]);      % Turning points
087 RFC = tp2rfc(TP);              % Rainflow cycles
088 paramD = [1 n n];
089 FrfcObs0 = cc2cmat(paramD,RFC); % Observed rainflow matrix
090 
091 FrfcObs = dtp2rfm(xD,n); % Observed rainflow matrix
092 
093 cmatplot(u,u,{FrfcObs/(T/2) Grfc},1)
094 
095 h=0.8; FrfcSmooth = smoothcmat(FrfcObs,1,h);
096 cmatplot(u,u,{FrfcObs/(T/2) FrfcSmooth/(T/2) Grfc},1)
097 
098 % Level Crossings
099 
100 mu = cmat2lc(param,Grfc);
101 muSum = cmat2lc(param,GrfcSum);
102 muObs = cmat2lc(param,FrfcObs/(T/2));
103 subplot(2,1,1), plot(mu(:,1),mu(:,2),muSum(:,1),muSum(:,2),'--')
104 subplot(2,1,2), plot(mu(:,1),mu(:,2),muObs(:,1),muObs(:,2),'--')
105 
106 % Damage
107 
108 beta = 4;
109 Dam = cmat2dam(param,Grfc,beta)
110 DamSum = cmat2dam(param,GrfcSum,beta)
111 DamObs = cc2dam(u(RFC),beta)/(T/2)
112 
113 Dmat = cmat2dmat(param,Grfc,beta);
114 DmatSum = cmat2dmat(param,GrfcSum,beta);
115 subplot(1,2,1), cmatplot(u,u,Dmat), axis('square'), v=axis;
116 subplot(1,2,2), cmatplot(u,u,DmatSum), axis('square'), axis(v)
117 
118 % Decomposition of a Mixed rainflow matrix
119 
120 % Model and Estimation
121 
122 n = 8; param = [-1 1 n];
123 M1.x0=[-0.4 -0.3]; M1.s=0.15; M1.lam=1; 
124 M2.x0=[0.3 0.4]; M2.s=0.15; M2.lam=1;
125 G1 = mktestmat(param,M1.x0,M1.s,M1.lam);
126 G2 = mktestmat(param,M2.x0,M2.s,M2.lam);
127 P=[1-0.1 0.1; 0.05 1-0.05]                % Transition matrix
128 [xD,z] = smctpsim(P,{G1 []; G2 []},5000); % Simulate
129 Fobs = dtp2rfm(xD,n);                     % observed mixed rainflow matrix
130 
131 known1.F = {G1 []; G2 []};   % known min-max and max-min matrices
132 init1.P = P;                 % initial guess of P-matrix
133 [Fest1,Est1] = estsmctp(Fobs,'P','ML',known1,[],init1);
134 Est1.P     % Estimated P-matrix
135 
136 known3.Ffun = 'f_funm';      % Function for calculating a submodel
137 known3.trModel2X = 'tr_m2x'; % transform from Model to X-vector
138 known3.trX2Model = 'tr_x2m'; % transform from X-vector to model
139 known3.param =param;
140 init3.P = P;       % initial guess of P-matrix
141 init3.M = {M1 M2}; % initial guess of Models for min-max mat
142 [Fest3,Est3] = estsmctp(Fobs,'P,CalcF','ML',known3,[],init3);
143 Est3.P     % Estimated P-matrix
144 Est3.M{:}  % Estimated parameters in models
145 
146 Pest_z = estmc(z,2)
147 
148 mc2stat(Est1.P)
149 mc2stat(Est3.P)
150 mc2stat(Pest_z)
151 mc2stat(P)
152 
153 % Methods for Estimation (Optional)
154 
155 known.F = {G1 []; G2 []};   % known min-max and max-min matrices
156 init.P = P;                 % initial guess of P-matrix
157 [Fest,Est] = estsmctp(Fobs,'P','ML',known,[],init); Est.P
158 [Fest,Est] = estsmctp(Fobs,'P','chi2',known,[],init); Est.P
159 [Fest,Est] = estsmctp(Fobs,'P','HD',known,[],init); Est.P
160 [Fest,Est] = estsmctp(Fobs,'P','KL',known,[],init); Est.P
161

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