Home > wafo > wdemos > itmkurs > itmkurs_lab4.m

# itmkurs_lab4

## PURPOSE

Script to computer exercises 4

## SYNOPSIS

This is a script file.

## DESCRIPTION

``` ITMKURS_LAB4 Script to computer exercises 4

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%```

## CROSS-REFERENCE INFORMATION

This function calls:
 cc2dam Calculates the total Palmgren-Miner damage of a cycle count. ccplot Plots a cycle count as a point process in the plane. cmat2dam Calculates the total Palmgren-Miner damage of a cycle matrix. cmat2dmat Computes the (Palmgren-Miner) damage matrix from a cycle matrix. cmatplot Plots a cycle matrix, e.g. a rainflow matrix. cocc Plots cycles as points together with isolines of a cycle matrix. dat2dtp The sequence of discretized turning points from a signal. dat2tp Extracts turning points from data, dcc2cmat Calculates the cycle matrix for a discrete cycle count. dtp2rfm Calculates rainflow matrix from discrete turning points. 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. 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 splitload Split a switching load (e.g. switchingload.mat) tp2mm Calculates min2Max and Max2min cycles from a sequence of turning points tp2rfc Finds the rainflow cycles from the sequence of turning points.
This function is called by:

## SOURCE CODE

```001 %ITMKURS_LAB4 Script to computer exercises 4
002 %
003 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
004 % Analysis of Measured Loads
005 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
006
007 %Estimation from Time Signal
008
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];
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];
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```

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