Home > wafo > cycles > rfmextrapolate.m

rfmextrapolate

PURPOSE ^

Extrapolates a rainflow matrix.

SYNOPSIS ^

[Fest,Pout,Fextreme,Fsmooth,Fest0] = rfmextrapolate2(F,Pin,plotflag)

DESCRIPTION ^

 RFMEXTRAPOLATE Extrapolates a rainflow matrix.
 
  CALL: Fest = rfmextrapolate(F);
        [Fest,Pout] = rfmextrapolate(F,Pin,plotflag);
        [Fest,Pout,Fextreme,Fsmooth,Fest0] = rfmextrapolate(F,Pin,plotflag)
 
    F           = Observed rainflow matrix                        [n,n]
    Pin         = Input parameters. (optional)           [struct array]
     .method    = Method for extrapolating the LC.  (see cmat2extralc for options)
                  'gpd,ml'  : Generalized Pareto distribution. (default)
                  'exp,mld' : Exponential distribution.  (Linear in lin-log)
     .u_lev     = Lower and upper levels for extrapolation of levelcrossings. 
                  (manual choice of levels)          u_lev=[i_min i_max]
     .LCfrac    = Fraction of LC for choosing thresholds for extraoplation.
                  (automatic choice of u_lev)       (default 0.05)
     .h         = Bandwidth for smoothing, h=-1 gives automatic choice.
     .beta      = Damage exponent.
     .Lim       = Limits where to use extreme RFM.  (manual choice)
                  Settings:  Lim.range  Lim.min  Lim.max  (see cmatcombine)
     .LimRelDam = Threshold for relative damage. 
                  (automatic choice of Lim)         (default 0.95)
     .param     = Defines discretization.
     .PL        = Values of countour lines.
    plotflag    = 0: Don't plot diagnostic plots, (default)
                  1: Plot final result of estimated limiting RFM,
                  2: Plot results of each step in the extrapolation.
 
    Fest        = Extrapolated RFM.                               [n,n]
    Pout        = Output parameters.                     [struct array]
    Fextreme    = Extreme RFM.                                    [n,n]
    Fsmooth     = Smoothed RFM.  (Kernel smoothing)               [n,n]
    Fest0       = Extrapolated RFM, before applying limits.       [n,n]
    
  Extrapolates the level crossing spectrum for high and for low levels. 
  Computes the 'extreme RFM', an approximation of the RFM which 
  is good for large cycles (i.e. for low minima and high maxima). 
  Computes the 'smoothed RFM' by using kernel smoothing.
  The final estimate of the 'extrapolated RFM' (or 'limiting RFM')
  is a combination of Fextreme and Fsmooth.
 
  Example:
    [G,Gh] = mktestmat([-1 1 64],[-0.2 0.2], 0.15,1);
    xD = mctpsim({G Gh},2000);
    Frfc = dtp2rfm(xD,64,'CS');
    Fest = rfmextrapolate(Frfc,[],1);
    Grfc = mctp2rfm({G Gh});
    cmatplot({Frfc Fest; Grfc G},4)
 
  See also  cmat2extralc, lc2rfmextreme, smoothcmat, cmatcombine

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

0001 function [Fest,Pout,Fextreme,Fsmooth,Fest0] = rfmextrapolate2(F,Pin,plotflag)
0002 
0003 %RFMEXTRAPOLATE Extrapolates a rainflow matrix.
0004 %
0005 % CALL: Fest = rfmextrapolate(F);
0006 %       [Fest,Pout] = rfmextrapolate(F,Pin,plotflag);
0007 %       [Fest,Pout,Fextreme,Fsmooth,Fest0] = rfmextrapolate(F,Pin,plotflag)
0008 %
0009 %   F           = Observed rainflow matrix                        [n,n]
0010 %   Pin         = Input parameters. (optional)           [struct array]
0011 %    .method    = Method for extrapolating the LC.  (see cmat2extralc for options)
0012 %                 'gpd,ml'  : Generalized Pareto distribution. (default)
0013 %                 'exp,mld' : Exponential distribution.  (Linear in lin-log)
0014 %    .u_lev     = Lower and upper levels for extrapolation of levelcrossings. 
0015 %                 (manual choice of levels)          u_lev=[i_min i_max]
0016 %    .LCfrac    = Fraction of LC for choosing thresholds for extraoplation.
0017 %                 (automatic choice of u_lev)       (default 0.05)
0018 %    .h         = Bandwidth for smoothing, h=-1 gives automatic choice.
0019 %    .beta      = Damage exponent.
0020 %    .Lim       = Limits where to use extreme RFM.  (manual choice)
0021 %                 Settings:  Lim.range  Lim.min  Lim.max  (see cmatcombine)
0022 %    .LimRelDam = Threshold for relative damage. 
0023 %                 (automatic choice of Lim)         (default 0.95)
0024 %    .param     = Defines discretization.
0025 %    .PL        = Values of countour lines.
0026 %   plotflag    = 0: Don't plot diagnostic plots, (default)
0027 %                 1: Plot final result of estimated limiting RFM,
0028 %                 2: Plot results of each step in the extrapolation.
0029 %
0030 %   Fest        = Extrapolated RFM.                               [n,n]
0031 %   Pout        = Output parameters.                     [struct array]
0032 %   Fextreme    = Extreme RFM.                                    [n,n]
0033 %   Fsmooth     = Smoothed RFM.  (Kernel smoothing)               [n,n]
0034 %   Fest0       = Extrapolated RFM, before applying limits.       [n,n]
0035 %   
0036 % Extrapolates the level crossing spectrum for high and for low levels. 
0037 % Computes the 'extreme RFM', an approximation of the RFM which 
0038 % is good for large cycles (i.e. for low minima and high maxima). 
0039 % Computes the 'smoothed RFM' by using kernel smoothing.
0040 % The final estimate of the 'extrapolated RFM' (or 'limiting RFM')
0041 % is a combination of Fextreme and Fsmooth.
0042 %
0043 % Example:
0044 %   [G,Gh] = mktestmat([-1 1 64],[-0.2 0.2], 0.15,1);
0045 %   xD = mctpsim({G Gh},2000);
0046 %   Frfc = dtp2rfm(xD,64,'CS');
0047 %   Fest = rfmextrapolate(Frfc,[],1);
0048 %   Grfc = mctp2rfm({G Gh});
0049 %   cmatplot({Frfc Fest; Grfc G},4)
0050 %
0051 % See also  cmat2extralc, lc2rfmextreme, smoothcmat, cmatcombine
0052 
0053 % References:
0054 %
0055 %   Johannesson, P., and Thomas, J-.J. (2001): 
0056 %   Extrapolation of Rainflow Matrices. 
0057 %   Extremes, Vol. 4, pp. 241-262.
0058 
0059 % Tested  on Matlab  5.3, 6.5
0060 %
0061 % History:
0062 % Created by PJ (Pär Johannesson) 24-Jul-2000
0063 % Updated by PJ 25-Apr-2003
0064 %   Added input parameter Pin.Lim
0065 
0066 % Check input arguments
0067 ni = nargin;
0068 no = nargout;
0069 error(nargchk(1,3,ni));
0070 
0071 if ni<2, Pin=[]; end
0072 if ni<3, plotflag=[]; end
0073 
0074 if isempty(plotflag)
0075   plotflag=0;
0076 end
0077 
0078 n = length(F);
0079 
0080 %
0081 % Default Parameters
0082 %
0083 
0084 %Pout.paramD = [1 n n];
0085 Pout.param = [1 n n];
0086 
0087 Pout.method = 'gpd,ml';
0088 Pout.u_lev = [];
0089 Pout.LCfrac = [];
0090 
0091 Pout.h = -1;   % Automatic choice of bandwidth for Kernel smoothing
0092 
0093 Pout.Lim = [];
0094 Pout.LimRelDam = [];
0095 Pout.beta = 7;
0096 
0097 % Plot-parameters
0098 Pout.PL = [10:20:90 99 99.9 99.99 99.999];
0099 
0100 % Copy input parameters Pin to Pout
0101 if ~isempty(Pin)
0102   Fname = fieldnames(Pin);
0103   for i = 1:length(Fname)
0104     Pout = setfield(Pout,Fname{i},getfield(Pin,Fname{i}));
0105   end
0106 end
0107 
0108 %
0109 % Default Parameters for extrapolatin of LC & where to use extreme RFM
0110 %
0111 
0112 if isempty(Pout.u_lev) & isempty(Pout.LCfrac)
0113     Pout.LCfrac = 0.05;
0114 end
0115 if isempty(Pout.Lim) & isempty(Pout.LimRelDam)
0116     Pout.LimRelDam = 0.95;
0117 end
0118 
0119 
0120 % 
0121 param = Pout.param;
0122 paramD = [1 n n];
0123 beta = Pout.beta;
0124 PL = Pout.PL;
0125 
0126 uD = levels(paramD);
0127 u = levels(param);
0128 du=u(2)-u(1);
0129 
0130 %
0131 % Extrapolate LC
0132 %
0133 
0134 % Find levels for extrapolation
0135 LCfrac = Pout.LCfrac;
0136 u_lev = Pout.u_lev;
0137 
0138 if plotflag >=2, figure, end
0139 
0140 if isempty(u_lev)
0141   slut=0;
0142   F1=F; i_lowPrev=0; i_highPrev=0;
0143   while ~slut
0144     lc = cmat2lc(paramD,F1);
0145     
0146     methodFindLev=3;
0147     if methodFindLev == 1
0148       [M,Mi] = max(lc(:,2)); Mi = Mi(1);
0149       dlc = [diff([0; lc(1:Mi-1,2)]); 0; flipud(diff([0; flipud(lc(Mi+1:end,2))]))];
0150       
0151       Nlc = sum(dlc~=0); % Number of non-zero in dlc = Number of jumps in lc
0152       Nextr = floor(0.2*Nlc); % Number of values in the tail of the distribution
0153       if Nextr<3, Nextr=3; end % At least 3 values in the tail of the distribution
0154       
0155       I = find(dlc~=0);
0156       i_low0 = I(Nextr+1);      % Low level
0157       i_high0 = I(end-Nextr);   % High level
0158       
0159     elseif methodFindLev == 2
0160       
0161       [M,Mi] = max(lc(:,2)); Mi = Mi(1);
0162       dlc1 = diff([0; lc(1:Mi-1,2)]);
0163       dlc2 = diff([0; flipud(lc(Mi+1:end,2))]);
0164       
0165       Nlc1 = sum(dlc1~=0); % Number of non-zero in dlc = Number of jumps in lc
0166       Nlc2 = sum(dlc2~=0); % Number of non-zero in dlc = Number of jumps in lc
0167       Nextr1 = max(ceil(0.4*Nlc1),3); % Number of values in the tail of the distribution
0168       Nextr2 = max(ceil(0.4*Nlc2),3); % Number of values in the tail of the distribution
0169       % At least 3 values in the tail of the distribution
0170       
0171       I1 = find(dlc1~=0);
0172       I2 = find(dlc2~=0);
0173       i_low0 = I1(Nextr1+1);      % Low level
0174       i_high0 = n-I2(Nextr2+1)+1;   % High level
0175       
0176       I = [I1; n-flipud(I2)+1];
0177       
0178     elseif methodFindLev == 3
0179       
0180       [M,Mi] = max(lc(:,2)); M = M(1); Mi=Mi(1);
0181       
0182       I = find(lc(:,2)>LCfrac*M);
0183       i_low0 = I(1);      % Low level
0184       i_high0 = I(end);   % High level
0185       
0186     end
0187     
0188     [M,Mi] = max(lc(:,2)); M = M(1); Mi=Mi(1);
0189     dlc = [diff([0; lc(1:Mi-1,2)]); 0; flipud(diff([0; flipud(lc(Mi+1:end,2))]))];
0190     Nlc = sum(dlc~=0); % Number of non-zero in dlc = Number of jumps in lc
0191     Nextr = 3; % Minimum number of values in the tail of the distribution
0192     I = find(dlc~=0);
0193     i_low = max(i_low0,I(Nextr+1));      % Low level
0194     i_high = min(i_high0,I(end-Nextr));   % High level
0195     
0196     % Remove cycles with  min>i_high  and  max<i_low
0197     F1 = F;
0198     for i = i_high+1:n
0199       F1(i,:) = 0;
0200     end
0201     for j = 1:i_low-1
0202       F1(:,j) = 0;
0203     end
0204     
0205     slut = (i_high==i_highPrev) & (i_low==i_lowPrev);
0206     i_highPrev=i_high; i_lowPrev=i_low;
0207     
0208   end
0209 else
0210   i_low = u_lev(1); i_high = u_lev(2);
0211   % Remove cycles with  min>i_high  and  max<i_low
0212   F1 = F;
0213   for i = i_high+1:n
0214     F1(i,:) = 0;
0215   end
0216   for j = 1:i_low-1
0217     F1(:,j) = 0;
0218   end
0219   lc = cmat2lc(paramD,F);
0220 end
0221 
0222 if plotflag >=2 % Diagnostic plot
0223   [M,Mi] = max(lc(:,2)); M = M(1); Mi=Mi(1);
0224   lc1 =lc; lc1(lc(:,2)==0,2) = 0.1;
0225   stairs(u(lc1(1:Mi,1)),lc1(1:Mi,2)), hold on
0226   stairs(u(lc1(Mi+1:end,1)-1),lc1(Mi+1:end,2)), hold on
0227   %stairs(lc(:,1),lc(:,2)+1), hold on
0228   %plot(lc(:,1),lc(:,2)+1), hold on
0229   %lc_min = min(lc(:,2));
0230   %I = find(dlc~=0);
0231   v = axis; axis([u([1 n]) 0.9 v(4)]);
0232   %plot(u(I),0.9*ones(1,length(I)),'*')
0233   plot(u([i_low i_high]),0.9*ones(1,2),'r.'), 
0234   plot(u([i_low i_low]),[0.9 v(4)],'r'),
0235   plot(u([i_high i_high]),[0.9 v(4)],'r'),
0236   if ~isempty(LCfrac), plot(u([1 n]),LCfrac*[M M],'r'), end
0237   hold off
0238   set(gca,'Yscale','log')
0239   title('Choice of threshold for extrapolation')
0240   xlabel('level'), ylabel('Number of crossings')
0241 end
0242 
0243 
0244 % Extrapolate LC
0245 
0246 u_lev = [i_low i_high];
0247 if plotflag>=2
0248   plotflag2 = 1; figure
0249 else
0250   plotflag2 = 0;
0251 end
0252 
0253 method = Pout.method;
0254 [lcEst,Est] = cmat2extralc(paramD,F,u_lev,method,plotflag2);
0255 
0256 if plotflag >= 2
0257   figure
0258   lcH =lcEst.High; lcH(lcH(:,2)==0,2) = 1e-100;
0259   lcL =lcEst.Low; lcL(lcL(:,2)==0,2) = 1e-100;
0260   plot(u(lc(:,1)),lc(:,2),'-'), hold on
0261   plot(u(lcH(:,1)),lcH(:,2),'r')
0262   plot(u(lcL(:,1)),lcL(:,2),'r') 
0263   hold off, grid on  
0264   set(gca,'YScale','log')
0265   Mlc = 10^(ceil(log10(max(lc(:,2))))); 
0266   v=axis; axis([u([1 n]) Mlc/1e8 Mlc])  
0267   title('Extrapolated level crossings')
0268   xlabel('level, u')
0269   ylabel('Level crossing intensity, \mu(u)')
0270 end
0271 Pout.u_lev = u_lev;
0272 
0273 %
0274 % Compute extreme RFM
0275 %
0276 
0277 lc1 = [0 0; lcEst.lc; n+1 0];
0278 
0279 Ga1 = lc2rfmextreme(lc1);
0280 Ga=Ga1(2:n+1,2:n+1);
0281 Ga(isnan(Ga)) = 0;
0282   
0283 if plotflag >= 2
0284   figure
0285   [qlGa PL] = qlevels(F,PL,u,u);
0286 
0287   contour(u,fliplr(u),flipud(F'),qlGa,'b'), hold on
0288   contour(u,fliplr(u),flipud(Ga'),qlGa,'r'), hold off
0289   title('Extreme RFM (red), compared with observed RFM (blue)')
0290   xlabel('min')
0291   ylabel('Max')
0292   v = axis; axis([min(u) max(u) min(u) max(u)]);
0293   axis('square')
0294 end
0295 
0296 %
0297 % Smooth the RFM
0298 %
0299 
0300 h=Pout.h;
0301 
0302 if h~=0
0303   % Find NOsubzero = Number of subdiagonals equal to zero
0304   i=0;
0305   while sum(diag(F,i)) == 0
0306     i=i+1;
0307   end
0308   NOsubzero = i-1;
0309   
0310   if h==-1 % Automatic choice of bandwidth
0311     h_norm = smoothcmat_hnorm(F,NOsubzero);
0312     % This choice is optimal if the sample is from a normal distribution
0313     % It usualy oversmooths, therefore one should choose a smaller bandwidth.
0314     
0315     h=0.5*h_norm; % Don't oversmooth !!!
0316   end
0317   
0318   [Gs] = smoothcmat(F,1,h,NOsubzero);
0319   
0320 else
0321   Gs=F;
0322 end
0323 
0324 Pout.h = h;
0325 
0326 if plotflag >= 2
0327   figure
0328   [qlGs PL] = qlevels(Gs,PL,u,u);
0329 
0330   contour(u,fliplr(u),flipud(Ga'),qlGs,'r'), hold on
0331   contour(u,fliplr(u),flipud(Gs'),qlGs,'b'), hold off
0332   title('Smooth RFM (blue), compared with Extreme RFM (red)')
0333   xlabel('min')
0334   ylabel('Max')
0335   v = axis; axis([min(u) max(u) min(u) max(u)]);
0336   axis('square')
0337 end
0338 
0339 % 
0340 % Combine 'Extreme RFM' and  RFMkernel
0341 %
0342 
0343 Dnorm = cmat2dam(paramD,F,beta)*100; % Suppose the measurement is 1/100 of the total lifetime
0344 ampF = cmat2amp(paramD,cmat2dmat(paramD,F/Dnorm,beta));
0345 ampGa = cmat2amp(paramD,cmat2dmat(paramD,Ga/Dnorm,beta));
0346 
0347 ptyp=2;
0348 Drel = zeros(n,1);
0349 IampF = (ampF(:,2)==0);
0350 IampGa = (ampGa(:,2)==0);
0351 I = ~IampF & ~IampGa;
0352 if ptyp < 3
0353   Drel(I)=ampGa(I,2)./ampF(I,2);
0354   Drel(IampF)=1e10;
0355 else
0356   Drel(I)=ampF(I,2)./ampGa(I,2);
0357   Drel(IampGa)=1e10;
0358 end
0359 
0360 Drel(IampF & IampGa)=0; 
0361 %Drel(isnan(Drel))=0; Drel(isinf(Drel))=1e10;
0362 
0363 if plotflag >= 2
0364   figure
0365   subplot(2,1,1),plot(ampF(:,1)*du,ampF(:,2)),hold on
0366   plot(ampGa(:,1)*du,ampGa(:,2),'r'),hold off
0367   v = axis; axis([0 n/2*du v(3:4)])
0368   title('Damage per amplitude'), xlabel('amplitude'), ylabel('D0, Dextreme')
0369   subplot(2,1,2),plot(ampF(:,1)*du,Drel), hold on
0370   if ptyp<3
0371     plot([0 n/2*du],0.95*[1 1],'r--'), hold off
0372     %v = axis; axis([0 n/2*du v(3:4)])
0373     v = axis; axis([0 n/2*du 0 5])
0374     title('Relative damage per amplitude'), xlabel('amplitude'), ylabel('Drel = Dextreme / D0')
0375     if ptyp==2, set(gca,'YDir','reverse'), end
0376   elseif ptyp==3
0377     plot([0 n/2*du],1/0.95*[1 1],'r--'), hold off
0378     v = axis; axis([0 n/2*du 0 5])
0379     title('Relative damage per amplitude'), xlabel('amplitude'), ylabel('Drel = D0 / Dextreme')
0380   end
0381   
0382 end
0383 
0384 Lim = Pout.Lim;
0385 
0386 if ~isfield(Lim,'range')
0387     % Where is the extreme RFM valid?
0388     % Choose limits according to damage.
0389     % At which amplitude does the extreme RFM give at least 95% of damage of F.
0390     LimRelDam = Pout.LimRelDam;
0391     I = find((ampF(:,2)~=0) & (ampGa(:,2)~=0));
0392     k = min(find(Drel(I)>LimRelDam));
0393     if ~isempty(k)
0394         Lim.range = I(k); % Valid fo amplitudes >= Lim.range
0395     else
0396         Lim.range=n; % Not valid, don't use Ga !!!
0397     end
0398 end
0399 
0400 [dummy,imax] = max(lcEst.lc(:,2));
0401 maxLC = imax(1);
0402 gap=floor(0.02*n)+1;
0403 if ~isfield(Lim,'min')
0404     Lim.min = maxLC-gap;
0405 end
0406 if ~isfield(Lim,'max')
0407     Lim.max = maxLC+gap;
0408 end
0409 
0410 [GaC0,LimOut] = cmatcombine(Ga,Gs,Lim);
0411 
0412 Pout.Lim = Lim;
0413 
0414 % Upper and lower bound of load values,
0415 % due to possible finite endpoints of the LC and hence of the RFM Ga
0416 
0417 f_max = sum(Ga);                   % distribution of maxima
0418 jj=min(find([0 fliplr(f_max)]>0)); % Find upper endpoint
0419 jmax = n-jj+2;                     % Largest possible maximum
0420 
0421 f_min = sum(Ga');          % distribution of minima
0422 ii=min(find([0 f_min]>0)); % Find lower endpoint
0423 imin = ii-1;               % Smallest possible minimum
0424 
0425 J=meshgrid(1:n);  % Columns = maximum
0426 I=J';             % Rows    = minimum
0427 K1 =  (I < imin) | (J > jmax); % Where shall the RFM be zero?
0428 GaC = GaC0;
0429 GaC(K1) = 0;
0430 
0431 Pout.imin = imin;
0432 Pout.jmax = jmax;
0433 
0434 if plotflag >= 1
0435   if plotflag>=2, figure, end
0436   [qlGaC PL] = qlevels(GaC,PL,u,u);
0437 
0438   contour(u,fliplr(u),flipud(F'),qlGaC,'b'), hold on
0439   contour(u,fliplr(u),flipud(GaC'),qlGaC,'r'), 
0440   %plot(u([1 n-LimOut.range]),u([LimOut.range n]),'g')
0441   plot(u([1 n-LimOut.range+1]),u([LimOut.range n]),'g')  % Mod PJ
0442   plot(u([1 n]),u([LimOut.max LimOut.max]),'g')
0443   plot(u([LimOut.min LimOut.min]),u([1 n]),'g'),hold off
0444   title('Limiting RFM (red), compared with observed RFM (blue)')
0445   xlabel('min')
0446   ylabel('Max')
0447   v = axis; axis([min(u) max(u) min(u) max(u)]);
0448   axis('square')
0449 end
0450 
0451 
0452 % Compare damage matrices
0453 if plotflag >= 2
0454   figure
0455   
0456   DmatGs = cmat2dmat(paramD,Gs,beta);
0457   DmatGa = cmat2dmat(paramD,Ga,beta);
0458   DmatGaC = cmat2dmat(paramD,GaC,beta);
0459   DmatF = cmat2dmat(paramD,F,beta);
0460   
0461   DamF = sum(sum(DmatF));
0462   DrelGa = sum(sum(DmatGa))/DamF;
0463   DrelGs = sum(sum(DmatGs))/DamF;
0464   DrelGaC = sum(sum(DmatGaC))/DamF;
0465 
0466   subplot(2,2,1), cmatplot(u,u,DmatF,13), title(['RFM, Drel = 1'])
0467   subplot(2,2,2), cmatplot(u,u,DmatGa,13), title(['Extreme RFM, Drel = ' num2str(DrelGa)])
0468   subplot(2,2,3), cmatplot(u,u,DmatGs,13), title(['Smoothed RFM, Drel = ' num2str(DrelGs)])
0469   subplot(2,2,4), cmatplot(u,u,DmatGaC,13), title(['Extrapolated RFM, Drel = ' num2str(DrelGaC)])
0470 
0471 end
0472 
0473 %
0474 % The Result!!!
0475 %
0476 
0477 Fest = GaC;
0478 
0479 if no>2
0480   Fest0 = GaC0;
0481   Fextreme = Ga;
0482   Fsmooth = Gs;
0483 end
0484 
0485

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