Post date: Mar 11, 2014 5:7:8 PM
///////////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
// Risk Calc version of the tank analysis
// See also the Risk Calc script "pressure.uc"
clear
_func beta() { // workaround until downs.cpp's beta can be fixed
if $1+$2+0 then a_ = a_; // checks dimensionlessness
if (($1===0) & ($2===0)) return([0,1]) // should have shape='beta'
if ($1===0) return(0) // shape='beta'
if ($2===0) return(1) // shape='beta'
if ($1 == 0) then if ($2 == 0) then return _env([0,1],_pbounds('beta',0,1,$1+1e-10,$2+1e-10)) else return _env(0,_pbounds('beta',0,1,$1+1e-10,$2+0))
if ($2 == 0) then return _env(1,_pbounds('beta',0,1,$1+0,$2+1e-10)) else _return _pbounds('beta',0,1,$1+0,$2+0);
}
func cb() return beta([$1, $1+1], [$2-$1, $2-$1+1])
s = cb(1,150)
ss = min(s,.1)
t = cb(0,2000)
tt = min(t,.005)
k2 = cb(3,500)
s1 = cb(0, 460)
k1 = cb(7,10000)
r = cb(0, 380)
e = t ||| (k2 ||| (s |&| (s1 ||| (k1 ||| r))))
ee = min(e1, 0.03)
// if failures of s1, k1 and r are dependent
e1 = t ||| (k2 ||| (s |&| (s1 | (k1 | r))))
ee1 = min(e1, 0.03)
// if failures of s, s1, k1 and r are dependent
ed = t ||| (k2 ||| (s & (s1 | (k1 | r))))
eed =min(ed, 0.03)
// if no assumption about independence can be made
ef = t | (k2 | (s & (s1 | (k1 | r))))
eef = min(ef, 0.03)
show ee; show ee1 in red; show eed in purple; show eef in blue
///////////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
// make up some interval values to play with
clear
function rand() {r__=random; return env(r__,random); }
a=rand(); b=rand(); c=rand(); d=rand(); e=rand()
//a;b;c;d;e
a = [ 0.03375244, 0.4723816]
b = [ 0.6900725, 0.9585024]
c = [ 0.7711181, 0.9777833]
d = [ 0.3543395, 0.6956177]
e = [ 0.09057617, 0.6003418]
////////////////////////////////////////////////////////////////////////////////
// "naive" calculation of the Birnbaum measure for the variables a
func F() return C - B * C * E - B * C * D + D * E - C * D * E - B * D * E + 2 * B * C * D * E
B = b
C = c
D = d
E = e
naive = F()
////////////////////////////////////////////////////////////////////////////////
// subinterval reconstitution is slow but rigorous
run brute
//Brute force interval library has been loaded
func map2() return C * (1 - B * ($1 + $2)) + $1 * $2 * (1 - C - B + 2 * B * C)
B = left(b)
C = right(c)
R1 = sir2(d,e,50)
B = right(b)
C = left(c)
R2 = sir2(d,e,50)
R = env(R1,R2)
R // [ 0.11755425, 0.7023882]
[ 0.1175542, 0.7023882]
////////////////////////////////////////////////////////////////////////////////
// Monte Carlo searching yields an inner estimate
function try() {
A = specific(a)
B = specific(b)
C = specific(c)
D = specific(d)
E = specific(e)
return C - B*C*E - B*C*D + D * E * (1 - C - B + 2 * B * C)
}
mc = try()
for i = 1 to 1000 do mc = env(mc, try())
////////////////////////////////////////////////////////////////////////////////
// algebraic rearrangements are rigorous, but limited
a0 = tri(c + d*e - b*c*e - c*d*e - b*d*e - b*c*d + 2*b*c*d*e)
a1 = tri(c - b*c*e - b*c*d + d * e * (1 -c - b + 2 * b * c))
a2 = tri(c * (1 - b * e - b*d) + d * e * (1 - c - b + 2 * b * c))
a3 = tri(c * (1 - b * (d + e)) + d * e * (1 - c - b + 2 * b * c))
a4 = tri(c - b*c*e + d * (e * (1 -c - b + 2 * b * c) - b*c))
a5 = tri(c - b*c*d + e * (d * (1 -c - b + 2 * b * c) - b*c))
a;b;c;d;e
[ 0.03375244, 0.4723816]
[ 0.6900725, 0.9585024]
[ 0.7711181, 0.9777833]
[ 0.3543395, 0.6956177]
[ 0.09057617, 0.6003418]
////////////////////////////////////////////////////////////////////////////////
// checking endpoints seems to work, but I don't quite understand why,
// nor can I believe it will work generally
func endpoint() if ($2==0) return left($1) else return right($1)
procedure assignall() begin
A = endpoint(a,$1)
B = endpoint(b,$2)
C = endpoint(c,$3)
D = endpoint(d,$4)
E = endpoint(e,$5)
end
// get an anchor for the env function
assignall(0,0,0,0,0)
f = F()
// search corners, but don't repeat for the variable a
i = 0 // for i=0 to 0 do
for j=0 to 1 do for k=0 to 1 do for l=0 to 1 do for m=0 to 1 do begin
assignall(i,j,k,l,m)
ff = F()
f = env(f,ff)
end
ends = f
////////////////////////////////////////////////////////////////////////////////
// combine endpoint checking with monotonicity
assignall(0,0,0,0,0)
f = F()
for n = 0 to 1 do begin
B = endpoint(b, n) // monotone increasing
C = endpoint(c, not n) // monotone decreasing
for l = 0 to 1 do for m = 0 to 1 do begin
D = endpoint(d,l)
E = endpoint(e,m)
ff = F()
f = env(f,ff)
end
end
monends = f
////////////////////////////////////////////////////////////////////////////////
// summarize and compare output intervals
show mc in red
show R
show a0; show a1 in olive; show a2 in blue; show a3 in green; show a4 in teal; show a5 in gray
ff = tri(monends)
show ff in green
range(monends)
[ 0.1258789, 0.699457]
range(ends)
[ 0.1258789, 0.699457]
R
[ 0.1175542, 0.7023882]
mc
[ 0.1605074, 0.6508642]
range(a5);range(a4);range(a3);range(a2);range(a1);range(a0)
[ -0.4162429, 1.059946]
[ -0.4354022, 1.149603]
[ -0.2326928, 1.267755]
[ -0.2326928, 1.267755]
[ -0.439358, 1.331206]
[ -1.185823, 1.894515]
naive
[ -1.185823, 1.894515]
///////////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////////
///////////////////////////////////////////////////////////////////////////////////////////
// Reliability of a noncoherent system
// Consider a pump, whose reliability is p, pumping liquid into
// a vessel monitored by a level sensor whose reliability is s,
// which informs a level controller whose reliability is c.
// The reliability of this system is nonunate. Below is an
// elaboration of the interval calculation of the system's
// reliability
//////////////////////////////////////////////////////////////////////////////////////
// naive calculation
func ncs() return $1 * $3 + $2 - $2 * $3
//////////////////////////////////////////////////////////////////////////////////////
// corners calculation
func lr() if ($2) return right($1) else return left($1)
func corners() begin
d = lr($3,0)
e = lr($1,0)*d+lr($2,0)*(1-d)
for i = 0 to 1 do for j = 0 to 1 do for k = 0 to 1 do begin
d = lr($3,k)
e = env(e, lr($1,i)*d+lr($2,j)*(1-d))
end
return e
end
//////////////////////////////////////////////////////////////////////////////////////
// subinterval reconstitution
run brute
//Brute force interval library has been loaded
func map3() return ncs($1,$2,$3)
clear
func sir() return sir3($1, $2, $3, $4)
//////////////////////////////////////////////////////////////////////////////////////
// Monte Carlo inner bounds
func mca() begin
mc_a = ncs(specific($1),specific($2),specific($3)) // this line should be replaceable by "mca = empty"
for mc_i = 1 to ($4-1) do begin
mc_a = env(mc_a,ncs(specific($1),specific($2),specific($3)))
// P = specific($1)
// S = specific($2)
// C = specific($3)
// mc_a = env(mc_a,mc__a)
// if (!contains(c,mc__a)) say P,S,C
end
return mc_a
end
//////////////////////////////////////////////////////////////////////////////////////
// translate input from experts into parameters
t = 500 hours
func pexp() return 1 - exp(-$2 * $1) // pexp(x, lambda) yield cumulative probability at x for exponential distribution with mean 1/lambda
func q() return pexp(t, $1+$2) / (1 + $2/$1) // q(lambda, mu) translates lambda and mu values ($1 and $2) to characterizations of the component unavailability
//////////////////////////////////////////////////////////////////////////////////////
// do the analysis
pump = q1 = q(1e-3 per hour, 0.1 per hour)
sensor = q2 = q(2e-3 per hour, 5e-2 per hour)
controller = q3 = q(3e-3 per hour, 1/60 hours) + 0
pump; sensor; controller
0.00990099
0.03846154
0.1525342
ncs(pump,sensor,controller)
0.03410508
//////////////////////////////////////////////////////////////////////////////////////
// characterize the uncertainty arbitrarily...plus or minus 30%
pump = measure(q1, 30%)
sensor = measure(q2, 30%)
controller = measure(q3, 30%)
range pump; range sensor; range controller
[ 0.003960396, 0.01584159]
[ 0.01538461, 0.06153847]
[ 0.06101367, 0.2440548]
n = ncs(pump,sensor,controller)
c = corners(pump,sensor,controller)
clear; s = sir(pump,sensor,controller, 20)
m = mca(pump,sensor,controller, 1000)
show s in blue; show n in black; show c in green; show m in red
range(n); range(s); range(c); range(m)
[ 0.0006075025, 0.06446601]
[ 0.01199703, 0.05903612]
[ 0.01259648, 0.05875033]
[ 0.01516103, 0.05562977]
//////////////////////////////////////////////////////////////////////////////////////
// natural (minimal) uncertainty, from significant digits
pump = q1 = q([1e-3 per hour], [0.1 per hour])
sensor = q2 = q([2e-3 per hour], [5e-2 per hour])
controller = q3 = q([3e-3 per hour], 1/[60 hours]) + 0
pump; sensor; controller
[ 0.003322259, 0.02912622]
[ 0.02654867, 0.05263158]
[ 0.1208633, 0.1853325]
[1e-3 per hour]; [0.1 per hour]
[ 0.0005, 0.0015] per hour
[ 0.05, 0.15] per hour
[2e-3 per hour]; [5e-2 per hour]
[ 0.0015, 0.0025] per hour
[ 0.045, 0.055] per hour
[3e-3 per hour]; [60 hours]
[ 0.0025, 0.0035] per hour
[ 55, 65] hours
n = ncs(pump,sensor,controller)
c = corners(pump,sensor,controller)
clear; s = sir(pump,sensor,controller, 20)
m = mca(pump,sensor,controller, 1000)
clear; show s in blue; show n in black; show c in green; show m in red
range(n); range(s); range(c); range(m)
[ 0.01719587, 0.05482085]
[ 0.02199165, 0.05004216]
[ 0.02224406, 0.04979065]
[ 0.0235281, 0.04926684]
//////////////////////////////////////////////////////////////////////////////////////
// the q() function is has repetitions; are we evaluating it correctly?
// pump
lambda = [1e-3 per hour]
mu = [0.1 per hour]
P = q(lambda, mu)
p1 = q(left(lambda),left(mu))
p2 = q(left(lambda),right(mu))
p3 = q(right(lambda),left(mu))
p4 = q(right(lambda),right(mu))
P; env(p1,p2,p3,p4)
[ 0.003322259, 0.02912622]
[ 0.003322259, 0.02912622]
// sensor
lambda = [2e-3 per hour]
mu = [5e-2 per hour]
P = q(lambda, mu)
p1 = q(left(lambda),left(mu))
p2 = q(left(lambda),right(mu))
p3 = q(right(lambda),left(mu))
p4 = q(right(lambda),right(mu))
P; env(p1,p2,p3,p4)
[ 0.02654867, 0.05263158]
[ 0.02654867, 0.05263158]
// controller
lambda = [3e-3 per hour]
mu = 1/[60 hour]
P = q(lambda, mu)
p1 = q(left(lambda),left(mu))
p2 = q(left(lambda),right(mu))
p3 = q(right(lambda),left(mu))
p4 = q(right(lambda),right(mu))
P; env(p1,p2,p3,p4)
[ 0.1208633, 0.1853325]
[ 0.1208752, 0.1853214]
//////////////////////////////////////////////////////////////////////////////////////
// the q() function seems monotone in each variable
// (even though it's not) because the exponential cdf
// in the function hardly affects the result at all
t = 500 hours
func pexp() return 1 - exp(-$2 * $1) // pexp(x, lambda) yield cumulative probability at x for exponential distribution with mean 1/lambda
func q() return pexp(t, $1+$2) / (1 + $2/$1) // q(lambda, mu) translates lambda and mu values ($1 and $2) to characterizations of the component unavailability
func qq() return 1 / (1 + $2/$1) // replace cdf value with "1"
pump = q(1e-3 per hour, 0.1 per hour)
sensor = q(2e-3 per hour, 5e-2 per hour)
controller = q(3e-3 per hour, 1/60 hours) + 0
pump;sensor;controller
0.00990099
0.03846154
0.1525342
pump = qq(1e-3 per hour, 0.1 per hour)
sensor = qq(2e-3 per hour, 5e-2 per hour)
controller = qq(3e-3 per hour, 1/60 hours) + 0
pump;sensor;controller
0.00990099
0.03846154
0.1525424
Pump = qq([1e-3 per hour], [0.1 per hour])
Sensor = qq([2e-3 per hour], [5e-2 per hour])
Controller = qq([3e-3 per hour], 1/[60 hours]) + 0
Pump;Sensor;Controller
[ 0.003322259, 0.02912622]
[ 0.02654867, 0.05263158]
[ 0.1208791, 0.1853361]
Pump = q([1e-3 per hour], [0.1 per hour])
Sensor = q([2e-3 per hour], [5e-2 per hour])
Controller = q([3e-3 per hour], 1/[60 hours]) + 0
Pump;Sensor;Controller
[ 0.003322259, 0.02912622]
[ 0.02654867, 0.05263158]
[ 0.1208633, 0.1853325]
+++
// this is a non-coherent system, so the lower bound is not obtained using the lower bounds of the components
// it works if both a and b are working, or if a is working and b isn't, or if both are not working
// and it doesn't work if neither of these three cases is true
a=0.7
b=0.75
a*b + a*(1-b) + (1-a)*(1-b)
0.775
a = [0.7, 0.8]
b = [0.75, 0.90]
a*b + a*(1-b) + (1-a)*(1-b)
[ 0.6149999, 0.9950001]
func logb( ) _return _ln($1)/_ln($2)
_func andc( ) _begin
s_=_tan(_pi*1 radian*(1-$3)/4)
_if (_abs(s_-1)<1e-8) _return ($1+0|0)|*|($2+0|0) _else _if (_abs(s_)<=1e-10) _return _min($1+0|0,$2+0|0) _else _if (_abs(s_+1)<1e-8) _return _max(0,$1+$2-1) _else _return logb(1+(_exp(($1+0|0)*_ln(s_))-1)|*|(_exp(($2+0|0)*_ln(s_))-1)/(s_-1),s_)
_end
_func orc( ) _begin
s_=_tan(_pi*1 radian*(1-$3)/4)
_if (abs(s_-1)<1e-8) _return 1-(1-($1+0|0))|*|(1-($2+0|0)) _else _if (abs(s_)<=1e-10) _return _max($1+0|0,$2+0|0) _else _if (abs(s_+1)<1e-8) _return _min(1,($1+0|0)+($2+0|0)) _else _return 1-logb(1+(_exp((1-($1+0|0))*_ln(s_))-1)|*|(_exp((1-($2+0|0))*_ln(s_))-1)/(s_-1),s_)
_end
_func orx() {s_ = 0; _for i_ = 1 to _args _do s_ = _min(1,s_ + _arg(i_)); _return s_; }
_func andx() _return 0
orx((a&b), (a & !b), (!a & !b))
[ 0.4499999, 1]
a & b
[ 0.4499999, 0.8]
a & !b
[ 0, 0.25]
a & (not b)
[ 0, 0.25]
!a & !b
[ 0, 0.25]
(not a) & (not b)
[ 0, 0.25]