20. Uncertainty and sensitivity#

20.1. Monte Carlo uncertainty analysis#

20.1.1. Model objects#

BioSTEAM streamlines uncertainty analysis with an object-oriented framework where a Model object samples from parameter distributions and reevaluates biorefinery metrics at each new condition. In essence, a Model object sets parameter values, simulates the biorefinery system, and evaluates metrics across an array of samples.


Model objects are able to cut down simulation time by sorting the samples to minimize perturbations to the system between simulations. This decreases the number of iterations required to solve recycle systems. The following examples show how Model objects can be used.

20.1.2. Create parameter distributions#

Let’s first learn how to create common parameter distributions using chaospy.

A triangular distribution is typically used when the parameter is uncertain within given limits, but is heuristically known to take a particular value. Create a triangular distribution:

from warnings import filterwarnings; filterwarnings('ignore')
from chaospy import distributions as shape
lower_bound = 0
most_probable = 0.5
upper_bound = 1
triang = shape.Triangle(lower_bound, most_probable, upper_bound)
Triangle(0, 0.5, 1)

A uniform distribution is used when the theoretical limits of the parameter is known, but no information is available to discern which values are more probable. Create a uniform distribution:

from chaospy import distributions as shape
lower_bound = 0
upper_bound = 1
unif = shape.Uniform(lower_bound, upper_bound)

A large set of distributions are available through chaospy, but generally triangular and uniform distributions are the most widely used to describe the uncertainty of parameters in Monte Carlo analyses.

20.1.3. Parameter objects#

Parameter objects are simply structures BioSTEAM uses to manage parameter values and distributions.

This section is just to get you familiar with Parameter objects. All the fields that a Parameter object can have are described below. Don’t worry if you don’t fully understand what each field does. The main idea is that we need to define the setter function that the Parameter object uses to set the parameter value to the element (e.g. unit operation, stream, etc.) it pertains to. We can also pass a distribution (i.e. a chaospy distribution) that will be accessible for Model objects to sample from. As for the name, units of measure, and the baseline value, these are all for bookkeeping purposes. BioSTEAM incorporates the name and units of measure when creating a DataFrame of Monte Carlo results and parameter distributions. Parameter objects are created by Model objects which implicitly pass both the system affected by the parameter, and the simulate function. But don’t worry about these last two fields, they are automatically added by the Model object when creating the parameter.

simulate: [function] Should simulate parameter effects.

system: [System] System associated to parameter.

name: [str] Name of parameter.

units: [str] Units of measure.

baseline: [float] Baseline value of parameter.

element: [object] Element associated to parameter.

setter: [function] Should set the parameter.

distribution: [chaospy.Dist] Parameter distribution.

Hopefully things will be become clearer as we start to create the parameter objects in the following sections…

20.1.4. Create a model object#

Model objects are used to evaluate metrics around multiple parameters of a system.

Create a Model object of the sugarcane biorefinery with internal rate of return and utility cost as metrics:

from biorefineries import sugarcane as sc
import biosteam as bst
solve_IRR = sc.tea.solve_IRR
total_utility_cost = lambda: sc.tea.utility_cost / 10**6 # In 10^6 USD/yr
metrics = (bst.Metric('Internal rate of return', sc.tea.solve_IRR, '%'),
           bst.Metric('Utility cost', total_utility_cost, '10^6 USD/yr'))
model = bst.Model(sc.sys, metrics)

The Model object begins with no parameters:

(No parameters)
metrics: Internal rate of return [%]
         Utility cost [10^6 USD/yr]

20.1.5. Add cost parameters#

A cost parameter changes cost requirements but not affect mass and energy balances.

Add number of fermentation reactors as a “cost” parameter:

R301 = bst.main_flowsheet.unit.R301 # The Fermentation Unit
    name='Number of reactors',
    element=R301, kind='cost',
    distribution=shape.Uniform(4, 10),
    hook=lambda N: int(round(N)) # Make sure value is an integer
def set_N_reactors(N):
    R301.N = N

The decorator uses the function to create a Parameter object and add it to the model:

parameters = model.get_parameters()
(<Parameter: [Fermentation-R301] Number of reactors>,)

Calling a Parameter object will update the parameter and results:

set_N_reactors_parameter = parameters[0]
print(f'Puchase cost at 5 reactors: ${R301.purchase_cost:,.0f}')
print(f'Puchase cost at 8 reactors: ${R301.purchase_cost:,.0f}')
Puchase cost at 5 reactors: $2,623,357
Puchase cost at 8 reactors: $3,137,039

Add the fermentation unit base cost as a “cost” parameter with a triangular distribution:

reactors_cost_coefficients = R301.cost_items['Reactors']
mid = reactors_cost_coefficients.n # Most probable at baseline value
lb = mid - 0.1 # Minimum
ub = mid + 0.1 # Maximum
@model.parameter(element=R301, kind='cost',
                 distribution=shape.Triangle(lb, mid, ub))
def set_exponential_cost_coefficient(exponential_cost_coefficient):
    reactors_cost_coefficients.n = exponential_cost_coefficient

Note that if the name was not defined, it defaults to the setter’s signature:

(<Parameter: [Fermentation-R301] Number of reactors>,
 <Parameter: [Fermentation-R301] Exponential cost coefficient>)

20.1.6. Add isolated parameters#

An isolated parameter does not affect Unit objects in any way.

Add feedstock price as a “isolated” parameter:

feedstock = sc.sugarcane # The feedstock stream
lb = feedstock.price * 0.9 # Minimum price
ub = feedstock.price * 1.1 # Maximum price
@model.parameter(element=feedstock, kind='isolated', units='USD/kg',
                 distribution=shape.Uniform(lb, ub))
def set_feed_price(feedstock_price):
    feedstock.price = feedstock_price

20.1.7. Add coupled parameters#

A coupled parameter affects mass and energy balances of the system.

Add fermentation efficiency as a “coupled” parameter:

@model.parameter(element=R301, kind='coupled',
                 distribution=shape.Triangle(0.85, 0.90, 0.95))
def set_fermentation_efficiency(efficiency):
    R301.efficiency = efficiency

Add crushing capacity as a “coupled” parameter:

lb = feedstock.F_mass * 0.9
ub = feedstock.F_mass * 1.1
@model.parameter(element=feedstock, kind='coupled',
                 distribution=shape.Uniform(lb, ub))
def set_crushing_capacity(capacity):
    feedstock.F_mass = capacity

20.1.8. Evaluate metric given a sample#

The model can be called to evaluate a sample of parameters.

Note that all parameters are stored in the model in the same order they were added:

parameters: Fermentation-R301 - Number of reactors
            Fermentation-R301 - Exponential cost coefficient
            Stream-Sugarcane - Feedstock price [USD/kg]
            Fermentation-R301 - Efficiency
            Stream-Sugarcane - Capacity
metrics: Internal rate of return [%]
         Utility cost [10^6 USD/yr]

Get dictionary that contain DataFrame objects of parameter distributions:

df_dct = model.get_distribution_summary()
Element Name Units Shape lower upper
0 Fermentation-R301 Number of reactors Uniform 4 10
1 Stream-Sugarcane Feedstock price USD/kg Uniform 0.0311 0.038
2 Stream-Sugarcane Capacity Uniform 3e+05 3.67e+05
Element Name Units Shape lower midpoint upper
0 Fermentation-R301 Exponential cost coefficient Triangle 0.4 0.5 0.6
1 Fermentation-R301 Efficiency Triangle 0.85 0.9 0.95

Evaluate sample:

model([8, 100000, 0.040, 0.85, feedstock.F_mass]) # Returns metrics (IRR and utility cost)
-  Internal rate of return [%]   0.0874
   Utility cost [10^6 USD/yr]     -18.2
dtype: float64

20.1.9. Monte Carlo#

Sample from a joint distribution, and simulate samples:

import numpy as np
N_samples = 200
rule = 'L' # For Latin-Hypercube sampling
np.random.seed(1234) # For consistent results
samples = model.sample(N_samples, rule)
    notify=50 # Also print elapsed time after 50 simulations
[50] Elapsed time: 9 sec
[100] Elapsed time: 17 sec
[150] Elapsed time: 26 sec
[200] Elapsed time: 35 sec

Although the system uses the last solution as an initial guess for the next, each scenario may be vastly different and there is no guarantee that this initial guess would be any better. Luckily, BioSTEAM can minimize perturbations to the system between simulations by pre-sorting the scenarios as well as use convergence data from past simulations to make better starting guesses for accelerated recycle convergence:

    optimize=True # Optimize simulation order
convergence_model = bst.ConvergenceModel(
    # Recycle loop prediction will be based on model parameters
    notify=50, # Also print elapsed time after 100 simulations
[50] Elapsed time: 9 sec
[100] Elapsed time: 17 sec
[150] Elapsed time: 25 sec
[200] Elapsed time: 32 sec

Notice how it took less time (~%5 less) to run the Monte Carlo simulations after optimizing the simulation order and/or predicting recycle loop convergence.

All data from simulation is stored in <Model>.table:

model.table # All evaluations are stored as a pandas DataFrame
Element Fermentation-R301 Stream-Sugarcane Fermentation-R301 Stream-Sugarcane -
Feature Number of reactors Exponential cost coefficient Feedstock price [USD/kg] Efficiency Capacity Internal rate of return [%] Utility cost [10^6 USD/yr]
0 4 0.491 0.0361 0.906 3.4e+05 0.13 -18.3
1 4 0.52 0.0316 0.922 3.07e+05 0.155 -16.4
2 5 0.55 0.0362 0.903 3.11e+05 0.121 -16.8
3 7 0.475 0.0375 0.908 3.64e+05 0.128 -19.6
4 7 0.506 0.0346 0.894 3.41e+05 0.136 -18.5
... ... ... ... ... ... ... ...
195 4 0.573 0.0341 0.9 3.22e+05 0.139 -17.3
196 5 0.465 0.0325 0.939 3.25e+05 0.16 -17.3
197 9 0.484 0.0358 0.896 3.14e+05 0.124 -16.9
198 8 0.459 0.0341 0.877 3.62e+05 0.136 -19.7
199 5 0.535 0.0334 0.886 3.35e+05 0.14 -18.1

200 rows × 7 columns

20.2. Sensitivity with Spearman’s rank order correlation#

Model objects also presents methods for sensitivity analysis such as Spearman’s correlation, a measure of monotonicity between variables:

df_rho, df_p = model.spearman_r()
print(df_rho['-', 'Internal rate of return [%]'])
Element            Parameter
Fermentation-R301  Number of reactors             -0.228
                   Exponential cost coefficient   0.0475
Stream-Sugarcane   Feedstock price [USD/kg]       -0.847
Fermentation-R301  Efficiency                      0.352
Stream-Sugarcane   Capacity                        0.338
Name: (-, Internal rate of return [%]), dtype: float64

Create a tornado plot of Spearman’s correlation between all parameters and IRR:

bst.plots.plot_spearman_1d(df_rho['-', 'Internal rate of return [%]'],
                           index=[i.describe() for i in model.parameters],
                           name='IRR [%]')
(<Figure size 640x480 with 3 Axes>,
 <Axes: xlabel="Spearman's correlation with IRR [%]">)

20.3. Single point sensitivity#

A quick way to evaluate sentivity is through single point sensitivity analysis, whereby a metric is evaluated at the baseline and at the lower and upper limits of each parameter. This method ignores the interactions between parameters and their distributions, but can help screen whether a system is sensitive to a given parameter. Model objects also facilitate this analysis:

baseline, lower, upper = model.single_point_sensitivity()
Element  Feature
-        Internal rate of return [%]   0.137
         Utility cost [10^6 USD/yr]    -17.9
dtype: float64

Element                                                                  -
Feature                                        Internal rate of return [%] Utility cost [10^6 USD/yr]
Element           Feature
Fermentation-R301 Number of reactors                                 0.138                      -17.9
                  Exponential cost coefficient                       0.135                      -17.9
Stream-Sugarcane  Feedstock price [USD/kg]                           0.157                      -17.9
Fermentation-R301 Efficiency                                          0.12                      -18.2
Stream-Sugarcane  Capacity                                           0.128                      -16.1

Element                                                                  -
Feature                                        Internal rate of return [%] Utility cost [10^6 USD/yr]
Element           Feature
Fermentation-R301 Number of reactors                                 0.135                      -17.9
                  Exponential cost coefficient                       0.138                      -17.9
Stream-Sugarcane  Feedstock price [USD/kg]                           0.115                      -17.9
Fermentation-R301 Efficiency                                         0.152                      -17.7
Stream-Sugarcane  Capacity                                           0.144                      -19.8

Create a tornado plot of the lower and upper values of the IRR:

IRR, utility_cost = model.metrics
metric_index = IRR.index
index = [i.describe(distribution=False) # Instead of displaying distribution, it displays lower, baseline, and upper values
         for i in model.parameters]
bst.plots.plot_single_point_sensitivity(100 * baseline[metric_index],
                                        100 * lower[metric_index],
                                        100 * upper[metric_index],
                                        name='IRR [%]',
(<Figure size 640x480 with 3 Axes>, <Axes: xlabel='IRR [%]'>)

Note that blue represents the upper limit while red the lower limit.