The Kitchin Research Group

Back in 2013 I wrote about using the uncertainties package to propagate uncertainties. The problem setup was for finding the uncertainty in the exit concentration from a CSTR when there are uncertainties in the other parameters. In this problem we were given this information about the parameters and their uncertainties.

The exit concentration is found by solving this equation:

\(0 = Fa0 - v0 * Ca - k * Ca**2 * V\)

So the question was what is Ca, and what is the uncertainty on it? Finding Ca is easy with fsolve.

from scipy.optimize import fsolve

Fa0 = 5.0
v0 = 10.0

V = 66000.0
k = 3.0

def func(Ca, v0, k, Fa0, V):
    "Mole balance for a CSTR. Solve this equation for func(Ca)=0"
    Fa = v0 * Ca     # exit molar flow of A
    ra = -k * Ca**2  # rate of reaction of A L/mol/h
    return Fa0 - Fa + V * ra

Ca, = fsolve(func, 0.1 * Fa0 / v0, args=(v0, k, Fa0, V))
Ca

0.0050000000000000001

The uncertainty on Ca is a little trickier. A simplified way to estimate it is:

\(\sigma_{Ca} = \sqrt{(dCa/dv0)^2 \sigma_{v0}^2 + (dCa/dv0)^2 \sigma_{v0}^2 + (dCa/dFa0)^2 \sigma_{Fa0}^2 + (dCa/dV)^2 \sigma_{V}^2}\)

We know the σ_i for each input, we just need those partial derivatives. However, we only have the implicit function we used to solve for Ca, and I do not want to do the algebra to solve for Ca. Luckily, we previously worked out how to get these derivatives from an implicit function using autograd. We just need to loop through the arguments, get the relevant derivatives, and accumulate the product of the squared derivatives and errors. Finally, take the square root of that sum.

import autograd.numpy as np
from autograd import grad

# these are the uncertainties on the inputs
s = [None, 0.1, 0.2, 0.05, 100]

S2 = 0.0

dfdCa = grad(func, 0)
for i in range(1, 5):
    dfdarg2 = grad(func, i)
    dCadarg2 = -dfdarg2(Ca, v0, k, Fa0, V) / dfdCa(Ca, v0, k, Fa0, V)
    S2 += dCadarg2**2 * s[i]**2

Ca_s = np.sqrt(S2)
print(f'Ca = {Ca:1.5f} +\- {Ca_s}')

Ca = 0.00500 +\- 0.00016776432898276802

That is the same uncertainty estimate the quantities package provided. One benefit here is I did not have to do the somewhat complicated wrapping procedure around fsolve that was required with uncertainties to get this. On the other hand, I did have to derive the formula and implement them. It worked fine here, since we have an implicit function and a way to get the required derivatives. It could take some work to do this with the exit concentration of a PFR, which requires an integrator. Maybe that differentiable integrator will come in handy again!

Copyright (C) 2018 by John Kitchin. See the License for information about copying.

org-mode source

Org-mode version = 9.1.14