[SOLVED] Brightway2: MonteCarloLCA with ecoinvent 3.10

With ecoinvent 3.10, the multiImpactMonteCarloLCA function posted by Chris Mutel and Pascal Lesage here produces Monte Carlo results that skew significantly higher than the LCA scores, and longer runs of 10,000 iterations sometimes include NaN results.

A while ago, I used the multiImpactMonteCarloLCA function with ecoinvent 3.8, and I received results that made perfect sense. However, more recently, using ecoinvent 3.10 (cutoff), I consistently receive Monte Carlo results that skew significantly higher than the LCA scores, and longer runs of 10,000 iterations sometimes include NaN results. I find the same thing happens with ecoinvent 3.91 (cutoff). I have held back from rolling back to ecoinvent 3.8 at this time. I am wondering if changes to the uncertainty metadata of ecoinvent 3.10 are incompatible with multiImpactMonteCarloLCA?

Below is some code to illustrate the problem, where I produce an LCA score and do a short run of 100 iterations for the GWP of 1 kg of chromium steel (to eliminate the possibility that the problem is with my inventory data). Thanks in advance for any diagnosis or workaround. (Although inconvenient, I could roll back my project to ecoinvent 3.8 if that is the most timely solution.)

ei310 = bw.Database('ecoinvent-3.10-cutoff')
chromiumSteelActivity = ei310.get('3a2caa153f4b39980c02f7a307874b83') # 'market for steel, chromium steel 18/8' (kilogram, GLO, None)
functional_unit = {chromiumSteelActivity: 1}
ipcc_2021_method = ('IPCC 2021', 'climate change', 'global warming potential (GWP100)')
method_list = [ipcc_2021_method]
lca = bw.LCA(functional_unit, ipcc_2021_method)
lca.lci()
lca.lcia()
print('GWP of 1 kg chromium steel =', lca.score)

# Code adapted from a Brightway2 seminar presented by Chris Mutel and Pascal Lesage:
# https://github.com/PoutineAndRosti/Brightway-Seminar-2017/blob/master/Day%201%20PM/Brightway%20and%20uncertainty.ipynb

def multiImpactMonteCarloLCA(functional_unit, method_list, iterations):
    # Create A and B matrices ready for Monte Carlo
    MC_lca = bw.MonteCarloLCA(functional_unit)
    MC_lca.lci()
    # Create dictionary for method name : characterization factor (C) matrices
    C_matrices = {}
    # Loop through method names and populate the C_matrices dictionary
    for method in method_list:
        MC_lca.switch_method(method)
        C_matrices[method] = MC_lca.characterization_matrix
    # Create results array with rows = number of methods, columns = number of iterations
    results = np.empty((len(method_list), iterations))
    # Loop through the Monte Carlo iterations and populate the results array
    for iteration in range(iterations):
        next(MC_lca)
        for method_index, method in enumerate(method_list):
            results[method_index, iteration] = (C_matrices[method]*MC_lca.inventory).sum()
    return results

# Confidence interval code by Xavier Guihot, posted on stackoverflow.com

def confidence_interval(data, confidence=0.95):
    dist = NormalDist.from_samples(data)
    z = NormalDist().inv_cdf((1 + confidence) / 2.)
    h = dist.stdev * z / ((len(data) - 1) ** .5)
    return dist.mean - h, dist.mean, dist.mean + h

# Find mean and 95% confidence interval for Monte Carlo results
test_results = multiImpactMonteCarloLCA(functional_unit, method_list, 100)
lower_95th, mean, upper_95th = confidence_interval(test_results[method_index])
print('Lower 95th =', lower_95th, 'Mean = ', mean, 'Upper 95th =', upper_95th)

Output from the above:
GWP of 1 kg chromium steel = 5.110700854620202
Lower 95th = 5.779195233025968 Mean =  5.967848380411753 Upper 95th = 6.156501527797538

Edit: I believe I have found a solution to this issue. Having imported ecoinvent 3.10 using the "newer" method (bw.import_ecoinvent_release("3.10", "cutoff")), I then left the ecoinvent-3.10-biosphere database alone, deleted ecoinvent-3.10-cutoff, and re-imported it using the "older" method below:

ei310 = bw.SingleOutputEcospold2Importer(fpei310, 'ecoinvent-3.10-cutoff', reparametrize_lognormals = True)
ei310.apply_strategies()
ei310.statistics()
ei310.write_database()

fpei10 is the local filepath to the datasets folder after downloading and unzipping ecoinvent 3.10_cutoff_ecoSpold02.7z manually. With this procedure, the "older" method recognises that a matching ecoinvent biosphere database is already in place, so doesn't try to reimport it, and the reparametrize_lognormals = True switch makes sure that the median value for activity impacts with lognormal uncertainty distributions matches the static impact results. Now I receive expected results. Ideally, the "newer" import method would have a reparametrize_lognormals option, which would avoid the need for the above workaround, and I will request this via the appropriate channel.

By the way, please accept my apology for naively posting a confidence interval function above that assumes the results follow a normal distribution!

I'll leave this query open for another day in case anyone has any additional or alternative information.