Mourad

Scott tries to evaluate Mourad's expressions using RAMAS Risk Calc.  Scott's entries into Risk Calc's Listener are shown here in unbolded black font.  The output from Risk Calc is shown in bold black font.  Interspersed commentary explaining the inputs and the outputs are shown in blue.

First, let's enter the basic definitions that Mourad provides.

beta_SU = beta

beta_SD = beta_S

Okay.  It looks fairly straightforward.  There might be some repeated variables to worry about; I can't tell yet because it depends on what parameters are uncertain and which aren't.

You didn't define beta_SU and beta_SD.

If your inputs are all precise probability distributions (or scalars), the result will be a precise probability distribution too, assuming all the variables are independent.

If some of PDFs and some are intervals (which is a special case of a fuzzy number), then you'll get a p-box as the result, right?

If some inputs are PDFs and some are fuzzy numbers, the result is what I call a hybrid number.  It is basically a nested stack of p-boxes.

Do you wanna go ahead and specify the inputs?  You can make them up.

I tried out the formula with made-up scalars.  Which should be distributions?  Which should be intervals (to start with)?

Cheers,

Scott

// inputs

N ~ K ~ numeric

T_1 = 1

MDT_sd = 1

MTTR_S = 1

MRT_S = 1

be = 0.3

be_S = 0.2

DC = 0.5

lambda_S = 1

be_SU = 0.4

be_SD = 0.6

// beta_SU and beta_SD were not defined

_func factorial() {F_ = 1; _for i_ = 1 _to $1 _do F_ = F_ * i_;  _return F_; };

function PFS() begin

  N = $1

  K = $2

  AKN = factorial(N)/factorial(N-K)

  lambda_SD = DC * lambda_S

  lambda_SU = (1 - DC) * lambda_S

  lambda_SDind = (1-be_S) * lambda_SD

  lambda_SUind = (1-be) * lambda_SU

  lambda_Sind = lambda_SDind + lambda_SUind

  one = AKN * lambda_Sind * MDT_sd

  two = (be_SU*lambda_SU + be_SD*lambda_SD) * MDT_sd

  P = 1

  for i = 1 to K-1 do P = P * (lambda_SUind * (T_1/(i+1) + MRT_S) + lambda_SDind * MTTR_S)

  return one * P + two

  end

// lambda_S is in both lambda_SU and lambda_SD

pfs = PFS(10,5)

pfs

    12486.85115

PFD Probability of failure on demand (unavailability of an item)

PFH Probability of failure per hour (unconditional failure intensity "failure frequency")

They are two metrics for the safety instrumented systems' performance in terms of Safety Integrity

PFS probability of failing safely

N = 3

K = 2

m = 6

DC = [0.4, 0.5, 0.6]

T_ST = uniform(1400 hour, 1500 hour)

lambda = lognormal2(-13.2, 0.43) * units('per hour')

theta = [0.4, 0.45, 0.55, 0.7]

beta_PT = [0.1, 0.15, 0.2]

beta_ST = [0.01, 0.07, 0.1]

Okay, here are the graphs of these various uncertain quantities that Risk Calc automatically draws.  (We omit the graphs for the scalar values, N, K and m.  They're just spikes at those values.)

Now let us compute the first derived variable.

lambda_DU = (1-DC) * lambda

Math Problem: incompatible operands

Oops!  It doesn't work.  What is the problem?  Let us look at the operands involved.

(1-DC)

    [ 0.4, 0.5, 0.6]

lambda

    ~lognormal(range=[6.11e-07,5.60e-06],  mean=2.02e-06,  var=8.36e-13) hour-1

One operand is a triangular fuzzy number and the other is a probability distribution.  In theory, those can be combined together to obtain a hybrid number, but these objects and this combination are not natively supported by Risk Calc.  However, Risk Calc has a library for doing the calculations.  So let's load that library and try to do the calculation.

run hybrid

Hybrid number arithmetic library has been loaded

H(1-DC, a)                   Put the first operand into the hybrid variable a

H(lambda, b)                 Put the second one into b

Hmultiply(a,b,c)            Multiply them together

Hdisplay

Hsay(c)

Best:  ~lognormal(range=[3.05e-07,2.80e-06],  mean=1.01e-06,  var=2.09e-13) hour-1

Base: ~(range=[2.4e-07,3.3e-06], mean=[0.8e-6,1.2e-6], var=[6.7e-14,5.3e-13]) hour-1

Hshow(c)

This hybrid number is displayed (with the Hshow procedure).  It consists of many nested p-boxes, only two of which appear in the graph (below) by default.  The tightest p-box at the top of the hybrid number is shown in gray, while the widest p-box at the base of the hybrid number is shown in black.

lambda_int_PT = (1-beta_PT) * (1 - theta) * lambda_DU

Math Problem: incompatible operands

Hmm.  Okay, it looks like there's going to be a lot of hybrid (i.e., fuzzy-probabilistic) operations here, and they are both more complex and more cumbersome to carry out in Risk Calc.  So, for now, let's pretend all of the fuzzy structures are really just intervals.  This will greatly simply the calculations and let us focus on other issues at the start.  We can return and redo the calculations with hybrid arithmetic once we're sure there are no other complicating issues.  So, first, let's replace the triangular and trapezoidal fuzzy numbers with their respective interval ranges.

DC = range DC

theta = range theta

beta_PT = range beta_PT

beta_ST = range beta_ST

Now we can do the basic calculations without worrying about hybrids.

lambda_DU = (1-DC) * lambda

As shown below, this results is just the same as the outer p-box at the base of the previously computed hybrid number.

When the inputs are just intervals and probability distributions, the results are p-boxes.  These are straightforwardly computed in Risk Calc and automatically graphed.

lambda_ind_PT = (1-beta_PT) * (1 - theta) * lambda_DU

lambda_ind_ST = (1-beta_ST) * theta * lambda_DU

lambda_ind_DU = lambda_int_PT + lambda_int_ST

lambda_CCF_PT = beta_PT * (1 - theta) * lambda_DU

lambda_CCF_ST = beta_ST * theta * lambda_DU

lambda_CCF_DU = lambda_CCF_PT + lambda_CCF_ST

These results including the graphs above may not be quite correct, because there are repeated variables in the definitions that we did not take account of.  But the calculations will at least be conservative bounds on the respective quantities.  Let's continue to see if we can compute (conservative) bounds on the risks of interest before we review the issue of repeated quantities and search for ways to ensure that the calculations are also best possible in addition to being rigorous.

First, we'll need a formula for the binomial coefficient function 'choose' which was not in the default  library for RAMAS Risk Calc, version 4.

func factorial() {F = 1; for i = 1 to $1 do F = F * i;  return F; };

function choose() return factorial($1)/(factorial($2)*factorial($1-$2));

A = choose(N, N-K+1)

The value of A is just 3.

for j = 0 to m-1 do S.j = A * (((lambda_ind_DU + (lambda_ind_PT * j)) * T_ST)^(N-K+2) - (lambda_ind_PT * j * T_ST)^(N-K+2))/(lambda_ind_DU * (N-K+2) * m * T_ST)

S = 0;  for j = 0 to m-1 do S = S + S.j

We can now compute PFD, the probability of failure on demand (unavailability of an item).

PFD = S + (T_ST/2)*(lambda_CCF_DU + (lambda_CCF_PT * (m-1)))

PFD

    ~(range=[-9.02e-05,0.0168],  mean=[0.00008,0.00272],  var=[0,0.000012])

//////////////////////////////////////////////////////////////////////////////////////////////////////////////////

//////////////////////////////////////////////////////////////////////////////////////////////////////////////////

//////////////////////////////////////////////////////////////////////////////////////////////////////////////////

// probability of failure per hour (unconditional failure intensity, "failure frequency")

B = N * choose(N-1,K-1) / (m * T_ST)

Disagreement between theoretical and observed variance

S = 0

for j = 0 to m-1 do S = S + (((lambda_ind_DU + (lambda_ind_PT * j)) * T_ST)^(N-K+1) - (lambda_ind_PT * j * T_ST)^(N-K+1))/(N-K+1)

PFH = B * S + lambda_CCF_DU

Math Problem: cannot handle distributions straddling zero

PFH

    PFH

//////////////////////////////////////////////////////////////////////////////////////////////////////////////////

//////////////////////////////////////////////////////////////////////////////////////////////////////////////////

//////////////////////////////////////////////////////////////////////////////////////////////////////////////////

// REVISITING probability of failure per hour (unconditional failure intensity, "failure frequency")

// symplify the notation a bit

d = lambda_ind_DU

p = lambda_ind_PT

t = T_ST

// then

// when N=3 and K = 2, the sum simplifies to

// sum = (dt + pjt)^2 - (pjt)^2 = md²t²+2dt²p(m-1)m/2

// which implies

PFH = lambda_CCF_DU + 3 * t * (d^2 + d * p * (m-1))

PFH

    ~(range=[8.58552e-09,8.85161e-07],  mean=[2.811584692e-08,0.000000285],  var=[0,5.5266665204895e-14]) hour-1

// this result may still be puffy or wrong, because lambda_CCF_DU has repetitions in theta, lambda,

// DC, etc., and d and p also have internal repetitions, some of which are shared with lambda_CCF_DU

// and each other

// where should we use independence assumptions?

// Which of the variables in your analysis can I assume are stochastically

// independent of one another?  Can I presume, for example, that T_ST

// and lambda are independent?  It is certainly okay (and in some ways

// easier and more convenient) to say that you don't know what the

// dependence is between these distributions.  Of course, all the variables

// you've specified cannot be independent of each other.  For instance,

// lambda_ind_PT and lambda_ind_ST both depend on the same variables,

// theta and lambda_DU, the latter of which at least has aleatory uncertainty.

// What do you wish to say about possible dependencies involving the fuzzy

// variables?

Scott asked Mourad:

Which of the variables in your analysis can I assume are stochastically independent of one another?

Can I presume, for example, that T_ST and lambda are independent?  It is certainly okay (and in some ways easier and more convenient) to say that you don't know what the dependence is between these distributions.

All variables you've specified cannot be independent of each other.  For instance, lambda_ind_PT and lambda_ind_ST both depend on the same variables, theta and lambda_DU, the latter of which at least has aleatory uncertainty.

What do you wish to say about possible dependencies involving the fuzzy variables?

Scott

Oh, and, also this question.  Are N and K always going to be 3 and 2 for your application?  Or are those just possible numbers that you're starting with?

Sc

I think there is a possibility to change N and K later, so it's better to keep using those general formulas but you can fix them at 3 and 2 if it helps you and simplifies the task.

All the variables in the left part of the table are independent, so for example T_ST and Lambda are totally independent and the same thing for all the other variables in that column. On the other side of the table, all the variables depend on the variables of the other side or even on each other.

I don't know if that's what you want to know, please let me know if not

By the way, hard luck for both of our soccer teams we got the same results :(

Regards,

Mourad