Use the calculated failure probability in the cost function (RBDO)

Soheil_Radfar · November 20, 2020, 9:38am

Hi everyone,
I use “Quantile Monte Carlo (QMC) with IP” to solve my RBDO problem.
I want to include the calculated probability of failure (Pf) at each optimization step in the cost function. In other words, the Pf value be a part of the cost function. How can this be done with UQLab commands?
I defined “myRBDO_QIP.Results.History.Constraints” as a separate variable, but that doesn’t work. Do you have a suggestion?

Here is my cost function:

function Y = uq_Ex02_RubBreakwat_cost(X)
% Definition of Val: Volume of Main Armor Layer

Val = 0.5*X(:,1).*(((X(:,1)+X(:,2)+X(:,3))./sin(X(:,9)))+...
    ((X(:,2)+X(:,3))./sin(X(:,9))))+(X(:,4).*X(:,1));

% Definition of Vul: Volume of Underlayer

Vfl = 0.5*X(:,2).*(((X(:,2)+X(:,3))./sin(X(:,9)))+...
    ((X(:,3))./sin(X(:,9)))) + X(:,5).*X(:,6) + X(:,2).*...
    (X(:,4)+(X(:,1)+X(:,2))/3)+X(:,2).*X(:,3)./sin(X(:,8));

% Definition of Vcl: Volume of Core Layer

Vcl = 0.5*X(:,3).*(2*X(:,7)+X(:,3).*cot(X(:,9))+X(:,3).*cot(X(:,8)))+...
    0.5*(4+X(:,6)+(X(:,2)./sin(X(:,8)))+((X(:,1)+X(:,2))./sin(X(:,9))));

Pf = myRBDO_QIP.Results.History.Constraints.Pf;
% Penalty Function for Controling Overtopping Failure
failurep1 = Pf(:,1);
if failurep1 > 0.001
    Penalty1 = 1000000 * (failurep1 - 0.001); 
else
    Penalty1 = 0;
end

% Penalty Function for Controling Failure due to Armor Instability
failurep2 = Pf(:,2); 
if failurep2 > 0.001
    Penalty2 = 1000000 * (failurep2 - 0.001);
else
    Penalty2 = 0;
end

% Penalty Function for Controling Toe Erosion Failure
failurep3 = Pf(:,3); 
if failurep3 > 0.001
    Penalty3 = 100000 * (failurep3 - 0.001);
else
    Penalty3 = 0;
end

% Calculation of Construction Cost
Y = 822*Val + 166*Vfl + 131*Vcl + Penalty1 + Penalty2 + Penalty3;

end

Thanks for your consideration

moustapha · November 22, 2020, 3:09pm

Hi @Soheil_Radfar,

Unfortunately the RBDO module was not designed to solve the case you’re presenting. You would have to make quite a few modifications in the module to get to use the Pf in the cost function.

If your computational model is not too expensive (i.e. you are not trying to use a surrogate model), I would suggest directly computing the failure probability within the cost function. You may have a look at the function uq_evalPfMC.m for an idea of how you could compute a failure probability corresponding to a given design solution using Monte Carlo simulation.
If you are using a gradient-based solver, remember to reinitialize the random number generator seed at each iteration.

Let me know if you need further clarifications or if you would really like to go down the road of modifying the RBDO module, then I may give you clues.

Cheers,
Moustapha

Soheil_Radfar · November 23, 2020, 3:52pm

Thanks @moustapha

My computational model is not expensive and I am using Quantile Monte Carlo. I understood your proposed solution and tried to implement it, but I have not yet been able to compute the failure probability corresponding to each optimization step (given design solution).

Soheil_Radfar · November 26, 2020, 6:37pm

Dear @moustapha,

Unfortunately, despite many tries, I was not able to calculate the probability of failure using uq_evalPfMC in the cost function.
Can you please help me in this regard? How can I compute a failure probability corresponding to a given design solution?

Thanks in advance.

moustapha · December 1, 2020, 11:10am

Hi @Soheil_Radfar

I will try to provide you a short workflow on how you can compute the failure probability for a given design variable d. I made a few simplifications but this is just to show you the basic idea.

You first need to identify in your inputs which ones are design variables (denoted d below) and which one are environmental variables (denoted Z).
Normally you have already defined the environmental variables in the RBDO set-up file. So you can already generate the samples Z ( I assume here a MC of size 1e4 but you can update accordingly).
Below is a simple example of a function that given a set of design values, will compute corresponding failure probabilities using crude Monte Carlo. The function takes the following inputs:

d : set of design variables to be evaluated
Z : Sampled environmental variables
DesignOpts: A structure similar to the ones used to generate an input object in UQLab. It should contain the field .Marginals.Type where for each design variable you will give the probability distribution
std is a vector with the standard deviation associated to each of the design variables. This is a vector with the same size as the number of random variables. If the random variables are constant, you can set the corresponding value to 0.
idxDes is the index of the design variables in the final variables X used to evaluate the limit-state

idxEnv is the index of the environmental variables in the final variables X used to evaluate the limit-state

  function Pf = EvalPf(d,Z,DesignOpts,std, idxDes, idxEnv)

% Loop over all the design samples
for ii = 1:size(d,1)

  % Loop over each dimension of the design variable:
  for jj = 1 : size(d(ii,:),1)

      % I assume a fixed standard definition is given for each design
      % variable
      DesignOpts.Marginals(jj).Moments = [d(ii,jj), std] ;
  end

      % Build an Input object
      myInput = uq_createInput(XIOpts,'-private') ;

      % Sample the random  variables locally
      D = uq_getSample(myInput,1e4) ;

  % Now recombine D and Z according to how they enter in your final
  % objective function
  X = zeros(1e4,length(idxDes)+length(idxEnv)) ;
  X(:,idxDES) = D ;
  X(:,idxEnv) = Z ;

  % Evaluate the limit-state function
  Y(ii,:) = myConstraint(X) ;
  % Get the corresponding failure probability
  Pf(ii,:) = sum(Y(ii,:)<0)/N ;

end
end

Once you have this you can plug this into the objective function you have defined above and replace the line

Pf = myRBDO_QIP.Results.History.Constraints.Pf;

with

Pf = EvalPf(d,Z,DesignOpts,std, idxDes, idxEnv)

I am aware this may not work exactly as shown for your specific case, but the idea is to show you the workflow and the set of ingredients you would need. Feel free to then adapt and build your own evalPf function.

Best,
Moustapa
EvalPf.m (926 Bytes)

Soheil_Radfar · December 2, 2020, 9:41am

Dear @moustapha,

Thanks for your valuable feedback.
I modified EvalPf.m based on your suggestions as follows:

function Pf = EvalPf(d,Z,DesignOpts,std, idxDes, idxEnv)

% Loop over all the design samples
for ii = 1:size(d,1)
    
    % Loop over each dimension of the design variable:
    for jj = 1 : size(d(ii,:),1)
        % I assume a fixed standard definition is given for each design
        % variable

        %% 3.1 Design variables
        %
        DesignOpts.Marginals(1).Name = 'a';  % Armor-layer height - X1
        DesignOpts.Marginals(1).Type = 'Gaussian';

        DesignOpts.Marginals(2).Name = 'b';  % Underlayer height - X2
        DesignOpts.Marginals(2).Type = 'Gaussian';

        DesignOpts.Marginals(3).Name = 'c';  % Core height - X3
        DesignOpts.Marginals(3).Type = 'Gaussian';

        DesignOpts.Marginals(4).Name = 'Cw';  % Crest width -X4
        DesignOpts.Marginals(4).Type = 'Gaussian';

        DesignOpts.Marginals(5).Name = 't'; % Toe thickness - X5
        DesignOpts.Marginals(5).Type = 'constant';

        DesignOpts.Marginals(6).Name = 'r'; % Toe berm width -X6
        DesignOpts.Marginals(6).Type = 'constant';

        DesignOpts.Marginals(7).Name = 'Bc'; % Core width - X7
        DesignOpts.Marginals(7).Type = 'constant';

        DesignOpts.Marginals(8).Name = 'al';  % Leeward slope - X8
        DesignOpts.Marginals(8).Type = 'Gaussian';

        DesignOpts.Marginals(9).Name = 'as';  % Seaward slope - X9
        DesignOpts.Marginals(9).Type = 'Gaussian';
        
        DesignOpts.Marginals(jj).Moments = [d(ii,jj), std] ;
    end
    
    
    % Build an Input object
    myInput = uq_createInput(XIOpts,'-private') ;
    % Sample the random  variables locally
    D = uq_getSample(myInput,1e4) ;
    % Now recombine D and Z according to how they enter in your final
    % objective function
    X = zeros(N,length(idxDes)+length(idxEnv)) ;
    X(:,idxDES) = D ;
    X(:,idxEnv) = Z ;
    
    % Evaluate the limit-state function
    Y(ii,:) = myConstraint(X) ;
    % Get the corresponding failure probability
    Pf(ii,:) = sum(Y(ii,:)<0)/N ;
end
end

But I have not yet been able to get the input data from the main RBDO solution process and the EvalPf works separately. What is the solutions for this?
P.S.: I attached all files, here:
RBDO.zip (4.3 KB)

moustapha · December 17, 2020, 9:54am

Hi @Soheil_Radfar,

I have just modified your files. The trick was to use the .Parameters field in defining the cost function within uq_Ex02_RubBreakwat.m (4.8 KB) . Please have a look at the Model User Manual or this Example for more details on how to define parameters (other than the basic random variables) to a model in UQLab.

I also made a few moddifictions in EvalPf.m (1.8 KB) and uq_Ex02_RubBreakwat_cost.m (1.4 KB), I let you have a look. Everything runs but you probably need to calibrate your problem as of now, there is no feasible solution (and the third component of the cost function returns imaginary numbers…). Also note that I used a sample size of 1e4 for computing Pf using Monte Carlo. Your target Pf being at 1e-3, you probably need to increase the sample set size to reduce the variance of the Pf estimate or use a different method.

Cheers,
Moustapha

Soheil_Radfar · December 18, 2020, 10:53am

Great!
This solution worked well for my case.

Thanks a lot dear @moustapha for your detailed answer and support.

All the best,
Soheil