Issues with Nataf transform implementation
In this post we present some issues related to the Mapping method implemented in Effective Quadrature, the Nataf transform.
Recalling the post Sensitivity Analisys with Effective Quadratures(https://discourse.equadratures.org/t/sensitivity-analysis-with-effective-quadratures),
we want to construct a model f = \approx f(M, S, V_0, k, P_0, T_a, T_0) taking into account the correlation between the incertain input variables.
from equadratures import *
import numpy as np
import matplotlib.pyplot as plt
# inputs
order_parameters = 3
mass = Parameter(distribution='uniform',lower=30., upper=60., order=order_parameters)
area = Parameter(distribution='uniform',lower=0.005, upper=0.02, order=order_parameters)
volume = Parameter(distribution='uniform', lower=0.002,upper=0.01, order=order_parameters)
spring = Parameter(distribution='uniform', lower=1000, upper=5000., order=order_parameters)
pressure = Parameter(distribution='uniform', lower=90000., upper=110000., order=order_parameters)
ambtemp = Parameter(distribution='uniform',lower=290.,upper=296., order=order_parameters)
gastemp = Parameter(distribution='uniform', lower=340., upper=360., order=order_parameters)
parameters = [mass,area,volume,spring,pressure,ambtemp,gastemp]
def piston(x):
mass, area, volume, spring, pressure, ambtemp, gastemp = x[0],x[1],x[2],x[3],x[4],x[5],x[6]
A = pressure*area + 19.62*mass - (spring*volume)/(1./area)
V = (area/(2.*spring))*(np.sqrt(A**2 + 4*spring*pressure*volume*ambtemp/gastemp))
C = 2*np.pi*np.sqrt(mass/(spring + area**2*pressure*volume*ambtemp/(gastemp*V**2)))
return C
Lets assume a strong correlation coefficient \rho = 0.85 between the Initial Gas Volume V_0 and the Filling Gas Temperature T_0, as well as between Athmosferic Pressure P_{atm} and the Ambient Temperature T_{amb}.
A milder correlation between the remaing parameters is assumed, setting the relative off-diagonal elements to \rho =0.3.
#P_0, athomesferic pressure, and T_a, ambient temperature, are correlated:
r_45= 0.85
#V_0, initial gas volume, and T_0, filling gas temperature, are correlated:
r_26= 0.85
# the other parameters have a smaller correlation:
r_ij= 0.3
# correlation matrix:
# 7x7: the off-diagonal elements quantify
# the correlation between input parameters
R = np.ones([7,7])*r_ij
R[4,5]= R[5,4]= r_45
R[2,6]= R[6,2]= r_26
The resulting Correlation Matrix which quantify the linear relationship between the inputs is:
print ('Correlation Matrix R:')
print (R)
[[0.3 0.3 0.3 0.3 0.3 0.3 0.3 ]
[0.3 0.3 0.3 0.3 0.3 0.3 0.3 ]
[0.3 0.3 0.3 0.3 0.3 0.3 0.85]
[0.3 0.3 0.3 0.3 0.3 0.3 0.3 ]
[0.3 0.3 0.3 0.3 0.3 0.85 0.3 ]
[0.3 0.3 0.3 0.3 0.85 0.3 0.3 ]
[0.3 0.3 0.85 0.3 0.3 0.3 0.3 ]]
As reported in the procedure explained in the post (https://discourse.equadratures.org/t/sensitivity-analysis-with-effective-quadratures)
in the first step of the heuristics we construct the polynomial inside the Independent space:
mybasis = Basis('total-order')
mypoly = Poly(parameters, mybasis, method='least-squares', sampling_args={'mesh':'tensor-grid','subsampling-algorithm':'qr', 'sampling-ratio':1.0})
The next step consists in the evaluation of the model in the transformed coordinates: to this aim we firstly use the Correlation class as follows;
# correlation class:
corr_obj = Correlations(mypoly,R)
corr_obj.set_model(piston)
poly_corr = corr_obj.get_transformed_poly()
# statistics:
print ('Mean and Variance for the Correlated case:',poly_corr.get_mean_and_variance())
The implementation of the Nataf transform leads to the issue reported to the line below: the Newton optimization method is not able to reach the convergence while calculating the $\textit{intermediate correlation matrix } R'. [1]
`---------------------------------------------------------------------------
RuntimeError Traceback (most recent call last)
in ()
1 # correlation class:
----> 2 corr_obj = Correlations(mypoly,R)
3 corr_obj.set_model(piston)
4 poly_corr = corr_obj.get_transformed_poly()
5 # statistics:
/home/irene/Documenti/EQuad/Effective-Quadratures/equadratures/correlations.py in init(self, poly, correlation_matrix, verbose)
54
55 if (self.D[i].name!=‘custom’) or (self.D[j].name!=‘custom’):
—> 56 rho = optimize.newton(check_difference, self.R[i,j], maxiter=50)
57 else:
58 res = optimize.least_squares(check_difference, R[i,j], bounds=(-0.999,0.999), ftol=1.e-03)
/usr/lib/python2.7/site-packages/scipy/optimize/zeros.py in newton(func, x0, fprime, args, tol, maxiter, fprime2, x1, rtol, full_output, disp)
341 if disp:
342 msg = "Failed to converge after %d iterations, value is s" (itr + 1, p)
→ 343 raise RuntimeError(msg)
344
345 return _results_select(full_output, (p, funcalls, itr + 1, _ECONVERR))
RuntimeError: Failed to converge after 50 iterations, value is nan
`
References
[1] Melchers, R. E., (1945) Structural Reliability Analysis and Predictions. John Wiley and Sons, second edition.
