ITMKURS_LAB4 Script to computer exercises 4 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Analysis of Measured Loads %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Calculates the total Palmgren-Miner damage of a cycle count. | |
Plots a cycle count as a point process in the plane. | |
Calculates the total Palmgren-Miner damage of a cycle matrix. | |
Computes the (Palmgren-Miner) damage matrix from a cycle matrix. | |
Plots a cycle matrix, e.g. a rainflow matrix. | |
Plots cycles as points together with isolines of a cycle matrix. | |
The sequence of discretized turning points from a signal. | |
Extracts turning points from data, | |
Calculates the cycle matrix for a discrete cycle count. | |
Calculates rainflow matrix from discrete turning points. | |
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 SMCTP. | |
Simulates a switching Markov chain of turning points, | |
Smooth a cycle matrix using (adaptive) kernel smoothing | |
Split a switching load (e.g. switchingload.mat) | |
Calculates min2Max and Max2min cycles from a sequence of turning points | |
Finds the rainflow cycles from the sequence of turning points. |
001 %ITMKURS_LAB4 Script to computer exercises 4 002 % 003 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 004 % Analysis of Measured Loads 005 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 006 007 %Estimation from Time Signal 008 009 load switchingload 010 plot(x0(:,1),x0(:,2)) 011 012 n = 32; param = [-1 1.2 n]; % Define discretization 013 u = levels(param); % Discrete levels 014 delta = u(2)-u(1) % Discretization step 015 TP = dat2tp(x0,delta); % Get turning points and rainflow filter 016 017 TP0 = dat2tp(x0); 018 plot(x0(:,1),x0(:,2),TP0(:,1),TP0(:,2),TP(:,1),TP(:,2)) 019 length(x0), length(TP0), length(TP) 020 021 RFC0 = tp2rfc(TP0); subplot(1,2,1), ccplot(RFC0) 022 RFC = tp2rfc(TP); subplot(1,2,2), ccplot(RFC) 023 beta = 6; % Damage exponent 024 Dam0 = cc2dam(RFC0,beta) % Damage 025 Dam = cc2dam(RFC,beta) % Damage, after rainflow filter 026 Dam/Dam0 027 028 dtp = dat2dtp(param,TP,delta); % Discretized turning points 029 tp = [dtp(:,1) u(dtp(:,2))']; % Discretized turning points 030 T = length(dtp); % Number of turning points 031 clf, plot(x0(:,1),x0(:,2),TP(:,1),TP(:,2),tp(:,1),tp(:,2)) 032 v=axis; hold on, 033 plot([v(1:2)],[u(2:end)'-delta/2 u(2:end)'-delta/2],'k:') 034 hold off, axis([v(1:2) param(1:2)]) 035 036 rfc = tp2rfc(tp); % Get rainflow cycles 037 dam = cc2dam(rfc,beta) % Damage, after discretization & rainflow filter 038 dam/Dam0 039 040 tz = [2 2; 46.40 1; 114.5 2; 161 1; 225.1 2; 270 1; 337.5 2; ... 041 384.8 1; 433.2 2; 600 1]; 042 [xxd,xd,z] = splitload(dtp,tz); 043 plot(xxd{1}(:,1),xxd{1}(:,2)) 044 plot(xxd{2}(:,1),xxd{2}(:,2)) 045 hmmplot(xd(:,2),z,xd(:,1),[1 2],'','',1) 046 047 % Estimation of the Subloads 048 049 dtp1 = dat2tp(xxd{1}); 050 [mM1,Mm1] = tp2mm(dtp1); 051 F1 = dcc2cmat(mM1,n); 052 dtp2 = dat2tp(xxd{2}); 053 [mM2,Mm2] = tp2mm(dtp2); 054 F2 = dcc2cmat(mM2,n); 055 cmatplot(u,u,{F1 F2}) 056 057 [G1s,h1] = smoothcmat(F1); 058 G1 = smoothcmat(F1,1,1.0,0); 059 [G2s,h2] = smoothcmat(F2); 060 G2 = smoothcmat(F2,1,0.8,0); 061 cmatplot(u,u,{G1s G2s; G1 G2}) 062 cmatplot(u,u,{F1 F2; G1 G2}) 063 064 %Estimation of the Regime Process} 065 066 N1 = length(dtp1), N2 = length(dtp2) 067 N12 = 4; N21 = 4; 068 p1=N12/N1; p2=N21/N2; 069 P = [1-p1 p1; p2 1-p2] % P-matrix 070 statP = mc2stat(P) % Stationary distribution 071 072 GG = {G1 []; G2 []}; 073 [Grfc,mu_rfc] = smctp2rfm(P,GG); 074 cocc(param,RFC,Grfc) 075 Frfc = dtp2rfm(dtp(:,2),n); 076 cmatplot(u,u,{Frfc Grfc*T/2}) 077 078 beta = 6; 079 Dam0, Dam, dam 080 damG = cmat2dam(param,Grfc,beta)*T/2 081 damG/Dam0 082 beta = 4; 083 Dmat = cmat2dmat(param,Frfc,beta); 084 DmatG = cmat2dmat(param,Grfc,beta)*T/2; 085 cmatplot(u,u,{Dmat,DmatG},3) 086 087 [xsim,zsim] = smctpsim(P,GG,T); 088 figure(1), hmmplot(u(xd(:,2))',z,1:length(xd),[1 2],'','',1) 089 figure(2), hmmplot(u(xsim)',zsim,1:T,[1 2],'','',1) 090 091 % Decomposition of a Mixed Rainflow Matrix 092 093 n = 16; param = [-1 1.2 n]; % Define discretization 094 u = levels(param); % Discrete levels 095 delta = u(2)-u(1) % Discretization step 096 TP = dat2tp(x0,delta); % Get turning points and rainflow filter 097 098 dtp = dat2dtp(param,TP,delta); % Discretized turning points 099 tp = [dtp(:,1) u(dtp(:,2))']; % Discretized turning points 100 T = length(dtp); % Number of turning points 101 rfc = tp2rfc(tp); % Get rainflow cycles 102 beta = 6; 103 dam = cc2dam(rfc,beta) % Damage, after discretization & rainflow filter 104 dam/Dam0 105 106 Frfc = dtp2rfm(dtp(:,2),n); % Observed rainflow matrix 107 108 tz = [2 2; 46.40 1; 114.5 2; 161 1; 225.1 2; 270 1; 337.5 2; ... 109 384.8 1; 433.2 2; 600 1]; 110 [xxd,xd,z] = splitload(dtp,tz); 111 hmmplot(xd(:,2),z,xd(:,1),[1 2],'','',1) 112 dtp1 = dat2tp(xxd{1}); 113 [mM1,Mm1] = tp2mm(dtp1); 114 F1 = dcc2cmat(mM1,n); 115 G1 = smoothcmat(F1,1,0.8,0); 116 dtp2 = dat2tp(xxd{2}); 117 [mM2,Mm2] = tp2mm(dtp2); 118 F2 = dcc2cmat(mM2,n); 119 G2 = smoothcmat(F2,1,0.6,0); 120 121 known1.F = {G1 []; G2 []}; % known min-max and max-min matrices 122 init1.P = P; % initial guess of P-matrix 123 warning off % Don't display warnings 124 [Fest1,Est1] = estsmctp(Frfc,'P','ML',known1,[],init1); 125 Est1.P % Estimated P-matrix 126 mc2stat(Est1.P) % Estimated stationary distribution 127 128 known3.Ffun = 'f_funm'; % Function for calculating a submodel 129 known3.trModel2X = 'tr_m2x'; % transform from Model to X-vector 130 known3.trX2Model = 'tr_x2m'; % transform from X-vector to model 131 known3.param = param; 132 init3.P = P; % initial guess of P-matrix 133 M1.x0=[0.1 0.1]; M1.s=0.15; M1.lam=2; % submodel 1 134 M2.x0=[0.5 0.7]; M2.s=0.1; M2.lam=1; % submodel 2 135 init3.M = {M1 M2}; % initial guess of Models for min-max matrices 136 137 OPTIONS(2) = 1e-1; % the termination tolerance for x; 138 OPTIONS(3) = 1e-1; % the termination tolerance for F(x); 139 [Fest3,Est3] = estsmctp(Frfc,'P,CalcF','ML',known3,[],init3); 140 Est3.P % Estimated P-matrix 141 mc2stat(Est3.P) % Estimated stationary distribution 142 Est3.M{:} % Estimated parameters in models 143 144 beta = 3:0.2:8; 145 Dam0 = cmat2dam(param,Frfc,beta)/(T/2); % Damage from load signal 146 FrfcEst1 = smctp2rfm(Fest1.P,Fest1.F); 147 Dam1 = cmat2dam(param,FrfcEst1,beta); % Damage, scenario 1 148 FrfcEst3 = smctp2rfm(Fest3.P,Fest3.F); 149 Dam3 = cmat2dam(param,FrfcEst3,beta); % Damage, scenario 3 150 plot(beta,Dam0,'b',beta,Dam1,'r',beta,Dam3,'g') 151 plot(beta,Dam1./Dam0,'r',beta,Dam3./Dam0,'g') 152
Comments or corrections to the WAFO group