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:
 cc2cmat Calculates the cycle count matrix from a cycle count. cc2dam Calculates the total Palmgren-Miner damage of a cycle count. cmat2dam Calculates the total Palmgren-Miner damage of a cycle matrix. cmat2dmat Computes the (Palmgren-Miner) damage matrix from a cycle matrix. cmat2lc Calculates the level crossings from a cycle matrix. cmatplot Plots a cycle matrix, e.g. a rainflow matrix. dat2tp Extracts turning points from data, dtp2rfm Calculates rainflow matrix from discrete turning points. estmc estimates transition matrix P from a time series of a Markov chain. estsmctp Estimate SMCTP model from an observed rainflow matrix. hmmplot plots a Hidden Markov Model. levels Calculates discrete levels given the parameter matrix. mc2stat Calculates the stationary distribution for a Markov chain. mctp2rfm Calculates the rainflow matrix for a MCTP. mctpsim Simulates a Markov chain of turning points 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.
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