Posts in category calcul

Uncertainity propagation and reliability analysis calculation functions

I made general calculation scripts based on a general model of input data. The syntax of the input data has to follow the following model of script (Attachments: model.py) :

Model of script for input models

  1. The user has to define several model parameters :
  • The number of input random variables;
  • The number of ouput random variables;
  • Finite differente step factor: in Open TURNS 0.12.3 the finite difference step is equal for every variable. In order to take into account the unit of the random variables, I use a finite difference step based on the mean value of each random variable and I define the step as: dx = k_fd*mean(x).

Others parameters can be defined.

...
# --------------------------------------------------------------
# Model parameters
# --------------------------------------------------------------

#Number of input random variables (NbVA)
NbVA=4
#Number of output variables of the physical model (NbVS)
NbVS=1
#Finite differente step factor
k_fd=0.001
#...


#Other parameters
binf=0.0
#...
...
  1. The user has to define his physical model thanks to :
  • wrappers ;
  • OpenTURNSPythonFunction which can use wrappers.

The user has to name his function : PhysModel.

...
# --------------------------------------------------------------
# Physical model
# --------------------------------------------------------------

#Declaration of wrapper(s)
#PhysModel=NumericalMathFunction("my_wrapper")

#Use of OpenTURNSPythonFunction
class myfunction(OpenTURNSPythonFunction):

	def __init__(self):
		OpenTURNSPythonFunction.__init__(self, NbVA, NbVS)

	def f(self, X):
		d=5.0/384.0*X[3]*X[0]**4/(X[1]*X[2])
		return [d]

PhysModel = NumericalMathFunction(myfunction())
...
  1. The user has to define his probabilistic model following several recommandations :
  • each marginal has to be named thanks to the method setName();
  • I only use the independent copula.

The user has to name his composed distribution : ProbaModel.

...
# --------------------------------------------------------------
# Probabilistic model
# --------------------------------------------------------------

# 1. definition of the marginals (1D distribution)

L0=LogNormal(5000.0,50.0,binf,LogNormal.MUSIGMA)
L0.setName('L')
L1=LogNormal(30000.0,4500.0,binf,LogNormal.MUSIGMA)
L1.setName("E")
L2=LogNormal(1.0e9,1.0e8,binf,LogNormal.MUSIGMA)
L2.setName("I")
L3=LogNormal(10.0,3.0,binf,LogNormal.MUSIGMA)
L3.setName("p")
#...

# 2. definition of the distribution collection

CollectionD=DistributionCollection(NbVA)
CollectionD[0]=Distribution(L0)
CollectionD[1]=Distribution(L1)
CollectionD[2]=Distribution(L2)
CollectionD[3]=Distribution(L3)
#...

# 3. Copula definition

IC = IndependentCopula(NbVA)

# 4. Composed distribution definition

ProbaModel = ComposedDistribution(CollectionD,Copula(IC))

# 5. Mean point definition

mean=ProbaModel.getMean()
...
  1. The user can define the types of calculation performed at the beginning of the script (Calculation 1/2). The corresponding functions will run thanks to the second part of the calculation definition (Calculation 2/2). The user has to define the character string related to the numbers of the calculation.
...
# --------------------------------------------------------------
# Calculation 1/2
# --------------------------------------------------------------
# 0: Single evaluation of the physical model at mean point

# 11: Quadratic cumul analysis
# 12: Uncertainity propagation using Monte Carlo simulations

# 121: Basic statistics on Monte Carlo results

# in the future
# 13: Polynomial chaos expansion

# 21: Reliability analysis using FORM
# 22: Reliability analysis using SORM

# 31: Reliability analysis using Monte Carlo simulations
# 32: Reliability analysis using latin hypercube sampling
# 33: Reliability analysis using directional simulations

# 211: Reliability analysis using importance sampling


Calculation=["0","12","121"]

for i in range(len(Calculation)):
	if (Calculation[i]=="11"):
		from fc_CQ import Calcul_CQ
	if (Calculation[i]=="12"):
		from fc_MC import Calcul_MC

	if (Calculation[i]=="121"):
		from fc_BasicStat import Calcul_BasicStat

...

# --------------------------------------------------------------
# Calculation 2/2
# --------------------------------------------------------------

#Modification of the finite difference step value
step=NumericalPoint(NbVA)
for i in range(NbVA):
	step[i]=mean[i]*k_fd
myGradient = NonCenteredFiniteDifferenceGradient(step,PhysModel.getEvaluationImplementation())
myHessian = CenteredFiniteDifferenceHessian(step,PhysModel.getEvaluationImplementation())
PhysModel.setGradientImplementation(myGradient)
PhysModel.setHessianImplementation(myHessian)

for i in range(len(Calculation)):
	if (Calculation[i]=="0"):
		Mean_evaluation=PhysModel(mean)
		for j in range(NbVS):
			print ""
			print "Evaluation of the physical model at mean point - output",j+1,"/",NbVS,":", Mean_evaluation[j]
	if (Calculation[i]=="11"):
		resCQ=Calcul_CQ(PhysModel,ProbaModel)
	if (Calculation[i]=="12"):
		resMC=Calcul_MC(PhysModel,ProbaModel)

	if (Calculation[i]=="121"):
		resBasicStat=Calcul_BasicStat(resMC)
...

The calculation functions are independent files, except Calcul_BasicStat. Each file contains a python function based on Open TURNS. Each function can be divided into 2 parts :

  • the calculation definition and run;
  • the display of the results.

Function Quadratic cumul analysis

Attachements: fc_CQ.py

The function requires:

  • a NumericalMathFunction or physical model
  • a ComposedDistribution or probabilistic model

Note : In Open TURNS 0.12.3, the calculation of the importance factors is only available for scalar NumericalMathFunction.

...
# Quadratic cumul analysis
	RV_ProbabilisticModel = RandomVector(Distribution(ProbabilisticModel))
	Function = RandomVector(PhysicalModel,RV_ProbabilisticModel)
	QuadraticC = QuadraticCumul(Function)

	if (NbVS==1):
		IF=NumericalPoint(NbVA)
		for i in range(NbVA):
			IF[i] = QuadraticC.getImportanceFactors()[i]

	Mean1=QuadraticC.getMeanFirstOrder()
	Mean2=QuadraticC.getMeanSecondOrder()

	CovG=QuadraticC.getCovariance()

	stdev=NumericalPoint(NbVS)
	for i in range(NbVS):
		stdev[i] = sqrt(CovG[i,i])
...

Function Uncertainity propagation using Monte Carlo simulations

Attachements: fc_MC.py

The function requires:

  • a NumericalMathFunction or physical model
  • a ComposedDistribution or probabilistic model
...
# Monte Carlo
	In = ProbabilisticModel.getNumericalSample(10000)
	Out = PhysicalModel(In)

	Moyenne = Out.computeMean()
	standarddeviation=Out.computeStandardDeviationPerComponent()
	Skewness = Out.computeSkewnessPerComponent()
	Kurtosis = Out.computeKurtosisPerComponent()
...

Function Basic statistics

Attachements: fc_BasicStat.py

The function requires:

  • a NumericalSample (result of uncertainity propagation using Monte Carlo simulations.

This kind of calculation has to be performed just after an uncertainity propagation using Monte Carlo simulations. This function will be able to be improved.

Last modifications : 12/06/2009