Problem with FERUM in performing FORM reliability analysis

Hello

I need to perform Reliability Analysis using FERUM to obtain the Probability of Failure using FORM. My main code (transcript) is a procedure in which the input is let's say X. This main procedure uses equation solving techniques such as fzero and fsolve to obtain some parameters. When I run my main code, it works just fine by itself using deterministic input parameters (without using it as a function to link it with FERUM). I do not know WHY! when it comes to FERUM and I make one or some of the inputs as Random Variables the code doesn't work and shows:

Error using fzero

Function value at starting guess must be finite and real.

but the main code works with anyvalue of X, but when I make X as a random variable giving it any kind of distribution and run FERUM the error shows up after one iteration.

I would really appreciate any help. I would be more than glad to share any Matlab code.

[Walter Roberson]

Yes, the code would help for testing.

[Benjamin]

The main code is "Demand Splitting.m" that works seperately. To run this code seperately, just comment function lines (comment line #1, #381, and uncomment line #15).

but in order to perform the reliability analysis using FERUM with the random variable of X = Qb :

Run "inputfile.m"

Type in Matlab command window: ferum

choose option 10 (FORM)

Enter

Note: in the case FERUM wasn't used in the computer before, FERUM toolbox can easily be downloaded and put the folder next to the Matlab files, then set a path for FERUM by Matlab>Home>Set Path>Add with subfolders> Select FERUM Folder>Move to Bottom>Save

[Walter Roberson]

The ferum I found from the official site does not have an option 10, and FORM is option 3.

Reference to non-existent field 'ig_max'.

Error in form (line 28)

i_max = analysisopt.ig_max;

Error in ferum (line 102)

[formresults] = form(1,probdata,analysisopt,gfundata,femodel,randomfield);

[Walter Roberson]

I made some small modifications to DemandSplitting so that it stopped displaying the results of each fsolve(), and I added in cross-checks to be sure that each fsolve() produced a sane result. I then ran, and everything went to completion.

RESULTS FROM RUNNING FORM RELIABILITY ANALYSIS

Number of iterations:

100

Reliability index beta:

NaN

Failure probability:

NaN

Number of calls to the limit-state function:

299

..............................................................................................

The following parameters are now available in your current workspace:

formresults.iter = Number of iterations

formresults.beta = Reliability index beta from FORM analysis

formresults.pf1 = Failure probability pf1

formresults.dsptu = Design point u_star

formresults.alpha = Alpha vector

formresults.dsptx = Design point in original space

formresults.imptg = Importance vector gamma

formresults.gfcn = Recorded values of the limit-state function

formresults.stpsz = Recorded step size values

formresults.e1 = Recorded values of g-criterion

formresults.e2 = Recorded values og u-criterion

formresults.Recorded_u = Recorded values of u (optional)

formresults.Recorded_x = Recorded values og x (optional)

formresults.dbeta_dthetaf = Sensitivities of beta index w.r.t. to distribution parameters

formresults.dpf_dthetaf = Sensitivities of probability of failure w.r.t. to distribution parameters

formresults.dbeta_drho = Sensitivities of beta index w.r.t. to correlation coefficients

formresults.dpf_drho = Sensitivities of probability of failure w.r.t. to correlation coefficients

formresults.delta = "Normalized" sensitivities of beta index w.r.t. to mean

formresults.eta = "Normalized" sensitivities of beta index w.r.t. standard deviation

formresults.nfun = Number of calls to the limit-state function

..............................................................................................

Elapsed time is 750.015813 seconds.

[Benjamin]

I do appreciate your help. Yes I am actually looking for the value of Failure Probability. So, what modification you made, OR could you please attach the modified files?

Thank you.

[Walter Roberson]

The change I made to ferum.m is: change

if exist('analysisopt')

   if isfield(analysisopt,'echo_flag')
      echo_flag = analysisopt.echo_flag;
   else
      echo_flag = 1;
   end
end

to

if exist('analysisopt')

   if isfield(analysisopt,'echo_flag')
      echo_flag = analysisopt.echo_flag;
   else
      echo_flag = 1;
   end
else
    echo_flag = 1;
end

[Walter Roberson]

I have attached my modified version of DemandSplitting . I made no change to the algorithm: all I did was change a number of calls of the form

[x] = fsolve(func2,x0); % Call solver

into the form

[x, fval, errorflag] = fsolve(func2,x0, opts); % Call solver
assert(errorflag >= 0, 'bad func2');

where

opts = optimoptions('fsolve', 'display', 'none');

This does not affect the calculation at all: it just prevents it from spitting out endless messages saying that the optimization worked.

[Benjamin]

Thank you for all the modifications. However, in the DemandSplitting.m file that you modified, the random variable, Qb, should be commented (line 15). Since Qb is the random variable and it is inside a function it should have a distribution (its parameters are in the inputfile.m). So if you just uncomment it and don't let Qb get a fix value we still get the same error!!!!

Qb =

NaN

Error using fzero (line 328)

Function value at starting guess must be finite and real.

Is there any solution for this?

[Walter Roberson]

No, there is no solution for that.

Your inputfile_BJ has marg indicating starting point 10 and mean 10, with type 1 and transform type 3. That is % Normal marginal distribution with z(i,:) = ( x(i,:) - mean ) / stdv where x(i,:) is the starting point. With starting point and mean both being 10, that is going to give 0, so the result of x_to_u starts out as 0.

Then inside form,

[ G, grad_g ] = gfun(lsf,x,grad_flag,probdata,analysisopt,gfundata,femodel,randomfield);

where grad_flag is ffd which was extracted from analysis_opt . Inside gfun, case ffd, and evaluator flag set at 'basic',

[ g, dummy ] = gfunbasic(lsf,allx(:,i),'no ',probdata,analysisopt,gfundata);

The 'no ' there becomes grad_flag inside of gfunbasic .

gfunbasic hits the final "else" case, so it does

G = eval(expression);

which executes gfun_BJ(Qb) where Qb = 10 because that is the starting point from marg .

gfun_BJ executes DemandSplitting(Qb) . That executes the fzero() on fun1 with Qb 10, and that works this first time. After a bunch of calculations with that Qb = 10, you get out the result that MaxFr = 44 . Then back in gfun_BJ that is subtracted from 50 to get g = 6 . This is returned as the first output of gfunbasic

Back in gfun the section for ffd, the 10 is incremented to 10.02 and gfunbasic is called again, calling gfun_BJ again calling DemandSplitting(Qb) where Qb is now 10.02 . But it turns out that after all of those calculations with Qb = 10.02, that MaxF = 44 again, even though the input was slightly different. I did not attempt to follow the computations in DemandSplitting .

Thus in the gfun section for ffd, g_a_step_ahead becomes 6 as well as g having been 6. grad_g(j,i) = (g_a_step_ahead - g)/h; and since those are both 6, that gives you a grad_g of 0.

This grad_g is returned back to form which then does

   alpha = -grad_G / norm(grad_G);                                      % Alpha vector

but grad_G is 0, so alpha here at the form level becomes 0 / 0 = nan. That pollutes a bunch of calculations, but is not the main problem.

form continues on to

d = search_dir(G,grad_G,u)

grad_G is 0 as explained, and u is 0 as I explained early on, because data starting point 10 minus data mean 10 = 0 . search_dir does

alpha = -grad_G / norm(grad_G);

and with grad_G being 0, that makes alpha nan at this search_dir level, resulting in a search direction of nan. Then back in form,

u_new = u + step * d;

and with d being the nan returned from search_dir, that makes u_new nan.

At that point, the code then tries to process location nan, which leads to trying to fzero() on fun1 with Qb nan and that is what leads fzero to give the error message about the initial execution not being finite that you are reporting the error message for.

In short, the code assumes that the gradient is never 0, and that turns out to be a bad assumption. Even in cases where the gradient is not logically 0, the code should have known about floating point round-off that can lead to gradients evaluating numerically as 0 to floating point precision even if they would be distinctly non-zero if calculated with sufficient precision.

The DemandSplitting code does a number of calculations and does max() . If any of the calculations is insensitive to the exact value of the input, then you run the risk of returning the same MaxFr for two different inputs, potentially leading to a gradient that is logically 0.

You could try changing your starting point in marg in your inputfile_BJ so that you might possibly get a non-zero gradient, but that would not solve the problem that the code has faulty assumptions about gradient never becoming 0.

The initial part of DemandSplitting is well commented with remarks such as "%yield strength of corroded tensile rebar" but the calculation part of it is almost uncommented.

[Benjamin]

Thank you very much for you complete and detailed answer. I was able to solve the problem by after your explanations by changing the starting point and ffd.para parameters.

Answer 1

[Walter Roberson]

[Moving down to be an Answer]

No, there is no solution for that.

Your inputfile_BJ has marg indicating starting point 10 and mean 10, with type 1 and transform type 3. That is % Normal marginal distribution with z(i,:) = ( x(i,:) - mean ) / stdv where x(i,:) is the starting point. With starting point and mean both being 10, that is going to give 0, so the result of x_to_u starts out as 0.

Then inside form,

[ G, grad_g ] = gfun(lsf,x,grad_flag,probdata,analysisopt,gfundata,femodel,randomfield);

where grad_flag is ffd which was extracted from analysis_opt . Inside gfun, case ffd, and evaluator flag set at 'basic',

[ g, dummy ] = gfunbasic(lsf,allx(:,i),'no ',probdata,analysisopt,gfundata);

The 'no ' there becomes grad_flag inside of gfunbasic .

gfunbasic hits the final "else" case, so it does

G = eval(expression);

which executes gfun_BJ(Qb) where Qb = 10 because that is the starting point from marg .

gfun_BJ executes DemandSplitting(Qb) . That executes the fzero() on fun1 with Qb 10, and that works this first time. After a bunch of calculations with that Qb = 10, you get out the result that MaxFr = 44 . Then back in gfun_BJ that is subtracted from 50 to get g = 6 . This is returned as the first output of gfunbasic

Back in gfun the section for ffd, the 10 is incremented to 10.02 and gfunbasic is called again, calling gfun_BJ again calling DemandSplitting(Qb) where Qb is now 10.02 . But it turns out that after all of those calculations with Qb = 10.02, that MaxF = 44 again, even though the input was slightly different. I did not attempt to follow the computations in DemandSplitting .

Thus in the gfun section for ffd, g_a_step_ahead becomes 6 as well as g having been 6. grad_g(j,i) = (g_a_step_ahead - g)/h; and since those are both 6, that gives you a grad_g of 0.

This grad_g is returned back to form which then does

   alpha = -grad_G / norm(grad_G);                                      % Alpha vector

but grad_G is 0, so alpha here at the form level becomes 0 / 0 = nan. That pollutes a bunch of calculations, but is not the main problem.

form continues on to

d = search_dir(G,grad_G,u)

grad_G is 0 as explained, and u is 0 as I explained early on, because data starting point 10 minus data mean 10 = 0 . search_dir does

alpha = -grad_G / norm(grad_G);

and with grad_G being 0, that makes alpha nan at this search_dir level, resulting in a search direction of nan. Then back in form,

u_new = u + step * d;

and with d being the nan returned from search_dir, that makes u_new nan.

At that point, the code then tries to process location nan, which leads to trying to fzero() on fun1 with Qb nan and that is what leads fzero to give the error message about the initial execution not being finite that you are reporting the error message for.

In short, the code assumes that the gradient is never 0, and that turns out to be a bad assumption. Even in cases where the gradient is not logically 0, the code should have known about floating point round-off that can lead to gradients evaluating numerically as 0 to floating point precision even if they would be distinctly non-zero if calculated with sufficient precision.

The DemandSplitting code does a number of calculations and does max() . If any of the calculations is insensitive to the exact value of the input, then you run the risk of returning the same MaxFr for two different inputs, potentially leading to a gradient that is logically 0.

You could try changing your starting point in marg in your inputfile_BJ so that you might possibly get a non-zero gradient, but that would not solve the problem that the code has faulty assumptions about gradient never becoming 0.

The initial part of DemandSplitting is well commented with remarks such as "%yield strength of corroded tensile rebar" but the calculation part of it is almost uncommented.