Risk Analysis and Hydrologic Design Under Uncertainty

Jorge A. Ramirez

First Order Second Moment Composite Risk Analysis and Advanced First Order Second Moment Composite Risk Analysis

The failure of a system ( e.g. , a hydrologic system) occurs when the load on the system (L) exceeds the capacity (C) of the system. In general, both C and L are uncertain and their uncertainty is described in terms of probability density functions. Thus, for given design capacity and design load, there is always a non-zero probability that the load will exceed the capacity. This probability represents the risk of failure of the system, and it is usually denoted by R. As a function of the joint probability density function of (C,L), the risk of failure R can be obtained as,

> R:=Int(Int(f(c,l),l=c..infinity),c=-infinity..infinity);

If C and L are mutually independent, the joint density is equal to the product of the marginal densities. Thus, R is,

R:=Int(Int(f(c)*f(l),l=c..infinity),c=-infinity..infinity);

It is convenient to work with a Performance Variable , Z. Two very useful definitions of the Performance Variable are the Safety Margin:

Z:=C-L;

and the natural logarithm of the Safety Factor:

> Z:=ln(C/L);

Observe that with these definitions, the evaluation of the risk of failure, R reduces to evaluating the probability that Z < 0; that is,

R:=Prob(Z<0);

Clearly, in order to estimate R, it is then necessary to obtain the probability density function of the Performance Variable , Z. This can be done analytically for certain types of marginal distributions and certain simple types of Z functions, using the so-called Derived Distribution Approach.

For example, if C and L are Normal , then, Z = C - L is also Normal. Thus, the risk of failure is R = F (- b ) where F(.) is the cumulative standard normal distribution. The parameter b defined below is known as the Reliability Index .

R:=Phi(-beta);

> beta:=mu[z]/sigma[z];

> mu[z]:=mu[C]-mu[L];

> var[z]=var[C]+var[L];

>

As an additional example, consider the case of both C and L log-Normal , then, Z = ln(C) -ln( L) is again Normal and the above expression can be used.. Thus, the risk of failure is R = F (- b ) where F(.) is the cumulative standard normal distribution and where:

R:=Phi(-beta);

> beta:=mu[z]/sigma[z];

> mu[z]:=mu[ln(C)]-mu[ln(L)];

> var[z]=var[ln(C)]+var[ln(L)];

However, in most cases of multivariate Z, the derived distribution approach becomes rapidly intractable analytically. In addition, in may be that the marginal distributions are not known. Thus, one has to resort to approximate methods. These methods are illustrated below. They are First Order Second Moment Analysis, and Advanced First Order Second Moment Analysis. In these methods, the function Z is linearized about an appropriate point, and its first and second moments are obtained from the linearized function. FOSM analysis linearizes the function Z=G(x) about the mean of x.

Let Z = G(x 1,.... x m ) = C - L = G 1 (x 1 ,...x n ) - G 2 (x n+1 ,...x m ). Then,

>

> mu[z]:=G(mu[x[k]]);

> sigma[z]:=sqrt(Sum((Diff('G(x[k])',x[k])*sigma[x[k]])^2,k=1..m));

> beta:=mu[z]/sigma[z]; R:=Phi(-beta);

However, systems fail at points other that the mean. Consequently, it is desirable to carry out the First Order Taylor Series expansion about a point on the failure surface. After all, it is precisely this point at which an evaluation of risk is required. However, the difficulty stems from the fact that the failure surface is unknown, and consequently, the First Order Taylor Series expansion remains a function of an unknown failure point. The methodology known as Advanced First Order Analysis of Uncertainty is illustrated below for a simple design example.

The First Order Taylor series expansion of G(.) about the failure point (y 1 ,...y m ) is:

> Z:=G(y[k])+Sum(Diff('G(y[k])',x[k])*(x[k]-y[k]),k=1..m);

Thus, in general the mean of Z is:

> mu[z]:=G(y[k])+Sum(Diff('G(y[k])',x[k])*(mu[x[k]]-y[k]),k=1..m);

and the variance of Z is:

> sigma[z]:=sqrt(Sum((Diff('G(y[k])',x[k])*sigma[x[k]])^2,k=1..m));

Finally, under AFOAU (i.e., linearization about a point on the failure surface), the risk of failure is obtained as R = F (- b ) where F(.) is the cumulative standard normal distribution.

> beta:=mu[z]/sigma[z];R:=Phi(-beta);

In evaluating beta above, if the true failure point, y,is known then G(y)=0, because by definition, Z=0 on the failure surface. The problem in carrying out the AFOAU as indicated is that the failure point is unknown. A simple iterative procedure is implemented below.

1) Assume a trial value of the failure point.

2) Obtain the value of b using the equation above.

3) Test if trial value of failure point is true failure point, that is, test G(y) = 0.

4) If trail value is not true failure point, then obtain a new trial value of failure point as indicated below and iterate.

> alpha[k]:=Diff('G(y[k])',x[k])*sigma[x[k]]/sigma[z];

>

> y[k]=mu[x[k]]-beta*alpha[k]*sigma[x[k]];

Example:

Storm sewer risk is a function of uncertainties from hydrologic and hydraulic sources. One of these sources is the effect of construction errors on the sewer slope. You have been commissioned to perform a risk evaluation of a storm sewer system. The city engineer is concerned with the construction error in setting the correct slope specified in the hydraulic design.
In performing your analysis, use the Rational formula in order to compute the flow that must be evacuated; and Manning's equation to compute the storm sewer capacity. The risk evaluation should account for uncertainties in factors affecting both system loading and system capacity. These include pipe diameter, d; sewer slope,s; Manning's coefficient, n; runoff coefficient, C; area, A; and precipitation intensity, i.

In your analysis performboth FOA and AFOA of uncertainty. Evaluate the risk assuming that both the load and the capacity are normally distributed.

Sewer capacity can be computed using Manning equation and the System load (flow that must be evacuated) can be computed using the Rational Formula.

The worksheet below presents the solution of the risk estimation problem using the techniques illustrated above.

> restart;with(student);with(stats);

Define Capacity Function, Qc and Load Function, Qo:

Manning Equation: Qc in cfs, Diameter d in ft

> Qc:=(n,d,S)->0.463/n*d^(8/3)*S^(1/2);

Rational formula: Qo in cfs, rainfall intensity i in in/h, and contributing area A in Acres.

> Qo:=(C,i,A)->C*i*A;

> SF:=Qc(n,d,S)/Qo(C,i,A);

Define performance variable Z = G(x i ).

> G:=(n,d,S,C,i,A)->0.463/n*d^(8/3)*S^(1/2)-C*i*A;

>

Define the so-called sensitivities of the function G(.):

> Dn:=Diff(G(n,d,S,C,i,A),n);

> Dd:=Diff(G(n,d,S,C,i,A),d);

> DS:=Diff(G(n,d,S,C,i,A),S);

> DC:=Diff(G(n,d,S,C,i,A),C);

> Di:=Diff(G(n,d,S,C,i,A),i);

> DA:=Diff(G(n,d,S,C,i,A),A);

> Dn:=(n,d,S)->-.463*d^(8/3)*sqrt(S)/(n^2);

> Dd:=(n,d,S)->1.234666667*d^(5/3)*sqrt(S)/n;

> DS:=(n,d,S)->.2315000000*d^(8/3)/(n*sqrt(S));

> DC:=(C,i,A)->-i*A;

> Di:=(C,i,A)->-C*A;

> DA:=(C,i,A)->-C*i;

Define the means of each one of the independent variables defining the performance variable Z.

The design specifications for this problem are given in terms of the means of n, d, S, C, i, and A, as given below:

nm:=0.015;dm:=4.9;Sm:=0.004;Cm:=0.825;im:=4.613;Am:=10;

Define the variances of the independent variables.

Uncertainty in the design values is quantified in terms of the standard deviation of the design variables, as given below.

> sn:=0.0008295;sd:=dm->0.01*dm;sS:=0.0008;sC:=0.058575;si:=0.7992;sA:=0.5;

The following procedure performs FOA. Results are graphed showing the Risk of Failure R as a function of the safety factor SF, for each of 6 different values of the slope.

> risk:=proc() local dm, SF, mu, var1, var2, var, beta, riskfo: fopen(risk, WRITE, TEXT): for dm from 4.9 by -0.1 while dm >=0.1 do SF:=Qc(nm,dm,Sm)/Qo(Cm,im,Am): mu:=G(nm,dm,Sm,Cm,im,Am): var1:=(Dn(nm,dm,Sm)*sn)^2+(Dd(nm,dm,Sm)*sd(dm))^2+(DS(nm,dm,Sm)*sS)^2: var2:=(DC(Cm,im,Am)*sC)^2+(Di(Cm,im,Am)*si)^2+(DA(Cm,im,Am)*sA)^2: var:= var1 + var2: beta:=mu/sqrt(var): riskfo:=statevalf[cdf,normald[0,1]](-beta): fprintf(risk, `%g %g %g \n`, SF, riskfo, beta): od:fclose(risk): end;

> risk();

The following procedure performs AFOA. Results are graphed showing the Risk of Failure R as a function of the safety factor SF, for each of 6 different values of the slope.

> risk:=proc()

> local dm, SF, mu, var1, var2, var, betam, riskm,

Failure point vector, and derivatives evaluated at the mean and at failure point:

> nt, dt, St, Ct, it, At, Dnt, Ddt, DSt, DCt, Dit, DAt,

Parameter vector a i

> an, ad, aS, aC, ai, aA,

> kk, beta,riska:

> fopen(error, WRITE, TEXT): fopen(risk, WRITE, TEXT): fopen(riskn, WRITE, TEXT):

For a range of diameters from 4.9 to 0.1

> for dm from 4.9 by -0.1 while dm >=2.0 do

Compute Safety Factor and perform FOAU

> SF:=Qc(nm,dm,Sm)/Qo(Cm,im,Am):

> mu:=G(nm,dm,Sm,Cm,im,Am):

> var1:=(Dn(nm,dm,Sm)*sn)^2+(Dd(nm,dm,Sm)*sd(dm))^2+(DS(nm,dm,Sm)*sS)^2: var2:=(DC(Cm,im,Am)*sC)^2+(Di(Cm,im,Am)*si)^2+(DA(Cm,im,Am)*sA)^2:

> var:=var1 +var2:

> betam:=mu/sqrt(var):

> riskm:=statevalf[cdf,normald[0,1]](-betam):fprintf(risk, `%g %g %g \n`, SF, riskm, betam):

Perform AFOAU: Select an arbitrary trial value of a point on the failure surface.

> nt:=SF*nm: dt:=SF*dm: St:=SF*Sm: Ct:=SF*Cm: it:=SF*im: At:=SF*Am: kk:=G(nt,dt,St,Ct,it,At):

Iterate until a true failure point is found.

> while abs(kk) >=0.000001 do

Evaluate the mean of Z and the derivatives of G(x) at the current trial value of the failure point.

> Dnt:=Dn(nt,dt,St): Ddt:=Dd(nt,dt,St): DSt:=DS(nt,dt,St): DCt:=DC(Ct,it,At): Dit:=Di(Ct,it,At): DAt:=DA(Ct,it,At): mu:=G(nt,dt,St,Ct,it,At)+Dnt*(nm-nt)+Ddt*(dm-dt)+DSt*(Sm-St)+DCt*(Cm-Ct)+Dit*(im-it)+DAt*(Am-At):

Evaluate the variance of Z at the current value of the failure point

> var:=(Dnt*sn)^2+(Ddt*sd(dm))^2+(DSt*sS)^2+(DCt*sC)^2+(Dit*si)^2+(DAt*sA)^2:fprintf(error, `%g %g %g %g %g %g %g %g \n`,dm,nt,dt,St,Ct,it,At,kk): fprintf(error, `%g %g %g %g %g %g %g %g \n`,Dnt,Ddt,DSt,DCt,Dit,DAt,mu,var):

Evaluate value of reliability index at current value of failure point

> beta:=mu/sqrt(var):

Evaluate parameter vector a i

> an:=Dnt*sn/sqrt(var): ad:=Ddt*sd(dm)/sqrt(var): aS:=DSt*sS/sqrt(var): aC:=DCt*sC/sqrt(var): ai:=Dit*si/sqrt(var): aA:=DAt*sA/sqrt(var):

Obtain new trial value of failure point

> nt:=abs(nm-beta*an*sn): dt:=abs(dm-beta*ad*sd(dm)): St:=abs(Sm-beta*aS*sS): Ct:=abs(Cm-beta*aC*sC): it:=abs(im-beta*ai*si): At:=abs(Am-beta*aA*sA):

> kk:=G(nt,dt,St,Ct,it,At):

>

> od:

> riska:=statevalf[cdf,normald[0,1]](-beta):

> fprintf(riskn,`%g %g %g %g %g \n`, SF, riskm, betam, riska, beta):

> od:

> fclose(riskn): fclose(error):

> end;

>

> risk();

>

>