Thread Subject:

Monte Carlo Simulation, Problem Vectorizing Data

Im using the Monte Carlo simulation method to estimate the probability of failure of a corroding box girder with time and in turn the reliability. The objective is to determine the reliability-time profile for the structure. The script below determines the probability of failure and the reliability at time t=0. I want these results for (t=0,1,2,3...,100). The geometry of the section determines the strength namely, thickness and the plastic modulus, these parameters are slowly decreasing with time due to the effects of corrosion.

This equation relates thickness with time: T = t0 - (CR)*(x)
Where, T = Thickness, t0 = initial thickness, CR = Corrosion Rate and x = time in years.
This is the equation where I am having problems vectorizing the data.

Here is the script:
                            ............................................
l = 20;
b = 1;
h = 1;
nfailure = 0;
nsamples = 1e6;
cr = 40e-6 + 4e-6*randn(1,1e6); %Corrosion Rate
mf = 650e3; %Mean force (Uniformly Distributed)
vf = 65e3^2; %Variance of the force (Uniformly Distributed)
muf = log((mf^2)/sqrt(vf+mf^2)); %Mean force (Lognorm Distributed)
sigmaf = sqrt(log(vf/(mf^2)+1)); %Variance of the Force(Lognorm Distributed)
F = lognrnd(muf,sigmaf,1,1e6);
%Lines 13 through 18 deal with generating array of thicknesses
mt = 0.014;
vt = 1.96e-6;
mut = log((mt^2)/sqrt(vt+mt^2));
sigmat = sqrt(log(vt/(mt^2)+1));
t = lognrnd(mut,sigmat,1,1e6); %Array of initial thicknesses (lognormally distributed)
b1 = b - 2*t;
h1 = h - 2*t;
Z = (b*h.^2/4 - b1.*h1.^2/4); %Plastic Modulus
S = F*l/4; %Applied Bending Moment(LOAD)
R = 275e6*Z; %Resistant Bending Moment(Plastic BM)
SM = R - S; %Safety Margin
for i = 1 : 1e6;
if SM(i) < 0;
    nfailure = nfailure + 1;
end;
end;
Reliability = (mean(SM))/(std(SM)) %Reliability Index
fprintf('Refer to standard normal table')
pfailure = nfailure/nsamples
                                 .............................................

-How can I use this script to execute the year on year calculations to determine the probability of failure and the reliability taking into account the corrosion with time?
-What functions would I use?
-And how can I determine the distributions of the various parameters at intervals of one year?


Any help at all would be great!
Thanks inadvance.

Hi Roger,

SNIP

The following cycle
> for i = 1 : 1e6;
> if SM(i) < 0;
> nfailure = nfailure + 1;
> end;
> end;
may be vectorized by a single command:
nfailure = sum(SM<0);

I let a specialist answers the rest questions.
Mira