Convergence des moments statistiques de la fonction d'Ishigami en fonction de la taille du plan d'expériences avec le chaos polynomial

Testé avec la version 0.12.3 d'OpenTURNS

# Chargement des librairies
from openturns import *
from math import *
from openturns_viewer import ViewImage
 
# Initialisation de la graîne pour les simulations de Monte-Carlo
RandomGenerator().SetSeed(0)
 
# Définition de la fonction factorielle
def factorielle(x):
	if x==0:
		return 1
	else:
		return x*factorielle(x-1)
 
# Définition des constantes intervenant dans le fonction d'Ishigami
A = 7.0
B = 0.1 
 
#Modele physique : fonction d'Ishigami
class myfunction(OpenTURNSPythonFunction):
	def __init__(self):
		OpenTURNSPythonFunction.__init__(self, 3, 1)
	def f(self, X):
		H = sin(X[0]) + A*(sin(X[1])**2) + B*(X[2]**4)*sin(X[0])
		Y=[H]
		return Y
 
function=myfunction()
g = NumericalMathFunction(function)
 
# Nombre de variables aléatoires d'entrées
dim = 3
 
#Modele probababiliste
mean=NumericalPoint(3)
mean[0]=0.0
mean[1]=0.0
mean[2]=0.0
 
L1=Uniform(-pi,pi)
L2=Uniform(-pi,pi)
L3=Uniform(-pi,pi)
 
Collection=DistributionCollection(dim)
Collection[0]=Distribution(L1)
Collection[1]=Distribution(L2)
Collection[2]=Distribution(L3)
 
Copule=IndependentCopula(dim)
Modelproba=ComposedDistribution(Collection,Copula(Copule))
 
# Definition des polynômes unidimensionels du chaos polynomial
legendre = OrthogonalUniVariatePolynomialFamily(LegendreFactory())
polyColl = PolynomialFamilyCollection(dim)
polyColl[0]=legendre
polyColl[1]=legendre
polyColl[2]=legendre
 
# Affichage d'un polynôme de la base
enumFunc = EnumerateFunction(dim)
MultiIndice_1 = enumFunc(0)
print "Polynome 1 de la base : ", MultiIndice_1
 
productPoly = OrthogonalProductPolynomialFactory(polyColl,enumFunc)
basis = OrthogonalBasis(productPoly)
 
mesure = basis.getMeasure()
################################################################	
# Calcul du nombre de termes dans le dvpt sur CP
# On fixe le degré du chaos polynomial
p=5
P=factorielle(dim+p)/(factorielle(dim)*factorielle(p))
print "Taille du chaos polynomial pour p=3 : ", P
 
 
################################################################	
nb_calculs=range(0,3)
 
moyenne=[0]
variance=[0]
taille_plan=[0]
taille_base_finale=[0]
for i in nb_calculs:
	sampleSize = P+10*i
	taille_plan.append(sampleSize)
	print "Taille du plan d experiences : ", sampleSize
	regressionStrategy = RegressionStrategy(MonteCarloExperiment(sampleSize))
	#regressionStrategy = RegressionStrategy(LHSExperiment(sampleSize))		
 
	#Fixed Strategy
	basisSize = P
	#print "Taille de la base", basisSize
	fixedStrategy= FixedStrategy(basis,basisSize)
	#polynomialChaosAlgorithm = FunctionalChaosAlgorithm(g,Distribution(Modelproba),AdaptiveStrategy(fixedStrategy),ProjectionStrategy(regressionStrategy))
 
	#SequentialStrategy:	
	sequentialStrategy=SequentialStrategy(basis,basisSize)	
	#polynomialChaosAlgorithm = FunctionalChaosAlgorithm(g,Distribution(Modelproba),AdaptiveStrategy(sequentialStrategy),ProjectionStrategy(regressionStrategy))
 
	#CleaningStrategy:		
 
	maxIndex = 3*P #Numero du dernier polynome Psi que l on peut utiliser
	basisSize_cleaning = P #ici le basisSise_cleaning est la taille max de la base toleree
	significanceFactor = 1.0e-3	
	cleaningStrategy = CleaningStrategy(basis,maxIndex,basisSize_cleaning , significanceFactor,True)	
	polynomialChaosAlgorithm = FunctionalChaosAlgorithm(g,Distribution(Modelproba),AdaptiveStrategy(cleaningStrategy),ProjectionStrategy(regressionStrategy))
 
 
	polynomialChaosAlgorithm.run()
 
	#result	
	polynomialChaosResult = polynomialChaosAlgorithm.getResult()	
 
	#Get the polynomial chaos coefficients	
	coefficients = polynomialChaosResult.getCoefficients()	
 
	#Get the metamodel as a NumericalMathFunction	
	metaModel = polynomialChaosResult.getMetaModel()	
 
	#Get the basis		
	basisCollection = polynomialChaosResult.getBasis()	
	print "Taille de la base finale : ", polynomialChaosResult.getBasis().getSize()
 
	taille_base_finale.append(polynomialChaosResult.getBasis().getSize())
	#Get the distribution of variables U	
	mesure = polynomialChaosResult.getMeasure()	
 
	#Get the composedmodel which is the model of the reduced variables U	
	composedModel = polynomialChaosResult.getComposedModel	
 
 
	# Calcul des moments a partir des coefficients
	print "Moyenne a partir des coefficients : " , polynomialChaosResult.getCoefficients()[0]
	print "Erreur sur la moyenne : " , (polynomialChaosResult.getCoefficients()[0]-3.5)/3.5
	moyenne.append(polynomialChaosResult.getCoefficients()[0])
	# Calcul de la variance a partir des coefficients
	variance_coeffs=polynomialChaosResult.getCoefficients().norm2()-	polynomialChaosResult.getCoefficients()[0]**2
	print  "variance a partir des coefficients : " , variance_coeffs
	VarianceTheorique = (A**2)/8.0 + (B)*(pi**4)/5.0 + (B**2)*(pi**8)/18.0 +1.0/2.0
	print  "Erreur sur la variance : " , (variance_coeffs-VarianceTheorique)/VarianceTheorique,"\n"
	#print  "ecart-type a partir des coefficients : " , sqrt(variance_coeffs)
	variance.append(variance_coeffs)
 
fic_cv=open('cv_Ishigami_cleaning_5','w')
for i in nb_calculs:
	fic_cv.write(str(taille_plan[i+1])+' '+str(taille_base_finale[i+1])+' '+str(moyenne[i+1])+' '+str(variance[i+1])+'\n')
fic_cv.close()
 
 
print moyenne
print variance
# Variance theorique
VarianceTheorique = (A**2)/8.0 + (B)*(pi**4)/5.0 + (B**2)*(pi**8)/18.0 +1.0/2.0
print "variance theorique : ", VarianceTheorique
 
fic_coeff=open('coefficients_Ishigami_cleaning_5','w')
nb_coeffs =  polynomialChaosResult.getCoefficients().getDimension()
affichage_coefficients=range(0,nb_coeffs)
for i in affichage_coefficients:
	#print enumFunc(i), polynomialChaosResult.getCoefficients()[i]
	fic_coeff.write(str(enumFunc(i))+' '+str(polynomialChaosResult.getCoefficients()[i])+'\n')
fic_coeff.close()

Comments

1. souchaud@… -- 2011-04-22 15:21

For a polynomial chaos calculation using the fixed strategy, the number of polynoms can be calculated in a simple way by the method :

P=EnumerateFunction(NbVA).getStrateCumulatedCardinal(order)

where NbVA is the number of input random variable, order is the order of the development of the chaos.