ITMKURS_LAB2 Script to computer exercises 2 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% lab2.m Switching Markov Loads and Rainflow Analysis %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Calculates the cycle count matrix from a cycle count. | |
Calculates the total Palmgren-Miner damage of a cycle count. | |
Calculates the total Palmgren-Miner damage of a cycle matrix. | |
Computes the (Palmgren-Miner) damage matrix from a cycle matrix. | |
Calculates the level crossings from a cycle matrix. | |
Plots a cycle matrix, e.g. a rainflow matrix. | |
Extracts turning points from data, | |
Calculates rainflow matrix from discrete turning points. | |
estimates transition matrix P from a time series of a Markov chain. | |
Estimate SMCTP model from an observed rainflow matrix. | |
plots a Hidden Markov Model. | |
Calculates discrete levels given the parameter matrix. | |
Calculates the stationary distribution for a Markov chain. | |
Calculates the rainflow matrix for a MCTP. | |
Simulates a Markov chain of turning points | |
Makes test matrices for min-max (and max-min) matrices. | |
Calculates the rainflow matrix for a SMCTP. | |
Simulates a switching Markov chain of turning points, | |
Smooth a cycle matrix using (adaptive) kernel smoothing | |
Finds the rainflow cycles from the sequence of turning points. |
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
Comments or corrections to the WAFO group