Lab 09: Monte Carlo Simulation

Lab 09: Monte Carlo Simulation#

Examples (by instructor)#

In this lab, you’ll practice Monte Carlo simulation from Chapter 14: sampling from input distributions, propagating uncertainty through a model, and summarizing the output distribution to estimate failure rates and reliability.

Example 1: Estimating π with Monte Carlo#

Throw \(N\) random darts at a unit square. The fraction landing inside the quarter-circle estimates \(\pi/4\).

import numpy as np
import matplotlib.pyplot as plt
from scipy import stats

Example 2: Process Reliability — Propagating Input Uncertainty#

A reactor product must have purity ≥ 97.5%. Residence time \(\tau \sim \mathcal{N}(120, 12^2)\) s and temperature \(T \sim \mathcal{N}(400, 8^2)\) K vary each batch. The purity model is:

\[\text{Purity} = 100 \times \left(1 - 0.2\,e^{-k(T)\,\tau}\right), \qquad k(T) = e^{-1550/T}\]

Run 50,000 simulated batches and estimate the failure rate.

Example 3: Convergence — How Many Samples Are Enough?#

The MC estimate improves as \(N\) grows. Error shrinks as \(\sim 1/\sqrt{N}\). Track the running failure rate estimate as batches accumulate.

Warm-Up: Syntax Practice#

Short exercises to get comfortable with the MC building blocks before the main problems.

Exercise 1 — Sample from a distribution.

Using np.random.default_rng(seed=10), draw 1,000 samples from \(\mathcal{N}(\mu=500,\; \sigma=20)\). Print the sample mean and std. Then print the fraction of samples exceeding 530.

import numpy as np

rng = np.random.default_rng(seed=10)

samples = rng.normal(500, 20, 1000)

print(f"Sample mean : {np.mean(samples):.2f}")
print(f"Sample std  : {np.std(samples, ddof=1):.2f}")
print(f"Fraction > 530 : {np.mean(samples > 530):.4f}")

Sample mean : 498.62
Sample std  : 19.90
Fraction > 530 : 0.0550

Exercise 2 — Pass a sample through a model.

Feed concentration \(C_0 \sim \mathcal{N}(2.0,\; 0.1^2)\) mol/L. After a first-order reaction with fixed \(k=0.5\) min⁻¹ and \(t=3\) min, the outlet concentration is:

\[C = C_0 \cdot e^{-k\,t}\]

Using seed=5 and \(N=5{,}000\) samples, compute the mean and std of \(C\). Print what fraction of outlet concentrations fall below 0.25 mol/L.

import numpy as np

k, t = 0.5, 3.0

rng = np.random.default_rng(___)
C0  = rng.normal(___, ___, ___)

C_out = ___

print(f"Mean C_out : {np.mean(C_out):.4f} mol/L")
print(f"Std  C_out : {np.std(C_out, ddof=1):.4f} mol/L")
print(f"Fraction below 0.25 mol/L : {np.mean(C_out < 0.25):.4f}")

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
Cell In[3], line 5
      1 import numpy as np
      3 k, t = 0.5, 3.0
----> 5 rng = np.random.default_rng(___)
      6 C0  = rng.normal(___, ___, ___)
      8 C_out = ___

File numpy/random/_generator.pyx:4957, in numpy.random._generator.default_rng()

File _pcg64.pyx:123, in numpy.random._pcg64.PCG64.__init__()

File bit_generator.pyx:535, in numpy.random.bit_generator.BitGenerator.__init__()

File bit_generator.pyx:307, in numpy.random.bit_generator.SeedSequence.__init__()

TypeError: SeedSequence expects int or sequence of ints for entropy not 

Exercise 3 — Two uncertain inputs.

Both \(k \sim \mathcal{N}(0.5,\; 0.05^2)\) and \(t \sim \mathcal{N}(3.0,\; 0.2^2)\) are uncertain (with fixed \(C_0 = 2.0\) mol/L). Using seed=8 and \(N=10{,}000\), sample both, compute \(C = C_0\,e^{-k\,t}\), and print mean, std, and fraction below 0.25 mol/L. Compare your answer to Exercise 2 — does uncertainty in \(k\) and \(t\) increase or decrease the spread?

import numpy as np

C0 = 2.0

rng = np.random.default_rng(___)
k_s = rng.___(___, ___, ___)
t_s = rng.___(___, ___,  ___)

C_out2 = ___

print(f"Mean C_out : {np.mean(C_out2):.4f} mol/L")
print(f"Std  C_out : {np.std(C_out2, ddof=1):.4f} mol/L")
print(f"Fraction below 0.25 mol/L : {np.mean(C_out2 < 0.25):.4f}")

Exercise 4 — Summarize output distribution.

Using the C_out2 array from Exercise 3, print the 5th and 95th percentiles and plot a histogram with a vertical line at the 5th percentile.

import numpy as np
import matplotlib.pyplot as plt


# Compute percentiles
p5  = np.___(C_out2, ___)
p95 = np.___(C_out2, ___)

print(f"5th  percentile : {p5:.4f} mol/L")
print(f"95th percentile : {p95:.4f} mol/L")

fig, ax = plt.subplots(figsize=(7, 4))
ax.hist(C_out2, bins=50, color='steelblue', edgecolor='white', alpha=0.8, density=True)
ax.axvline(p5, color='tomato', linestyle='--', linewidth=2, label=f'5th pct = {p5:.3f}')
ax.set_xlabel('Outlet concentration (mol/L)')
ax.set_ylabel('Density')
ax.set_title('MC Output Distribution')
ax.legend()
plt.tight_layout()
plt.show()

Practice Problems (by students)#

Problem 1: Heat Exchanger — Outlet Temperature Reliability#

A shell-and-tube heat exchanger heats a process stream. The outlet temperature is modeled as:

\[T_\text{out} = T_\text{in} + \frac{U \cdot A \cdot \Delta T_{lm}}{\dot{m} \cdot C_p}\]

For simplicity, use the linearized form:

\[T_\text{out} = T_\text{in} + \eta \cdot (T_\text{hot} - T_\text{in})\]

where \(\eta\) is the heat exchanger effectiveness. The following inputs are uncertain each day:

Variable	Symbol	Distribution
Inlet temperature (°C)	\(T_\text{in}\)	\(\mathcal{N}(25,\; 3^2)\)
Hot side temperature (°C)	\(T_\text{hot}\)	\(\mathcal{N}(120,\; 5^2)\)
Effectiveness	\(\eta\)	\(\mathcal{N}(0.75,\; 0.04^2)\)

The process requires \(T_\text{out} \geq 80°C\).

(a) Using seed=15 and \(N = 50{,}000\) samples, simulate \(T_\text{out}\) for all three uncertain inputs. Print mean, std, 5th percentile, and the failure rate (% of days below 80°C).

import numpy as np
import matplotlib.pyplot as plt

# --- Model ---

# (a) All three inputs uncertain

(b) Plot the output histogram with vertical lines at the spec limit (80°C) and the mean.

# (b) Histogram with spec line and shaded failing region

(c) Which input contributes most to variability? Re-run the simulation three times, each time holding two inputs fixed at their means and only sampling the third. Print the std of \(T_\text{out}\) for each case and identify the most influential input.

# (c) One-at-a-time sensitivity

Problem 2: Batch Reactor Yield — Operating Condition Optimization#

A batch reactor produces a polymer. The yield (%) is modeled as:

\[Y = 100 \times \left(1 - e^{-k(T)\,\tau}\right), \qquad k(T) = A\,e^{-E_a/(R\,T)}\]

with \(A = 2\times10^6\) min⁻¹, \(E_a/R = 6{,}000\) K, and the following uncertain operating conditions:

Variable	Symbol	Baseline	Std dev
Residence time (min)	\(\tau\)	30	3
Temperature (K)	\(T\)	370	10

The minimum acceptable yield is 85%.

(a) Using seed=99 and \(N=50{,}000\), simulate the yield distribution at the baseline conditions. Print mean, std, and failure rate.

import numpy as np
import matplotlib.pyplot as plt

# --- Model ---

# (a) Baseline: tau~N(30,3), T~N(370,10)

(b) The process engineer proposes raising the mean temperature to 380 K (same std = 10 K). Re-run the simulation with this change. By how many percentage points does the failure rate drop?

# (b) Higher mean T: T~N(380,10)

(c) Alternatively, the engineer could tighten temperature control to \(\sigma_T = 5\) K (keeping mean at 370 K). Re-run this scenario. Compare the failure rate to part (b) — which strategy is more effective?

# (c) Tighter sigma_T: T~N(370,5)