* Using autograd for error propagation
:PROPERTIES:
:categories: autograd,uncertainty
:date: 2018/11/05 21:04:21
:updated: 2018/11/05 21:04:21
:org-url: http://kitchingroup.cheme.cmu.edu/org/2018/11/05/Using-autograd-for-error-propagation.org
:permalink: http://kitchingroup.cheme.cmu.edu/blog/2018/11/05/Using-autograd-for-error-propagation/index.html
:END:
Back in [[http://kitchingroup.cheme.cmu.edu/blog/2013/03/07/Another-approach-to-error-propagation/][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.
| Parameter | value | \sigma |
|-----------+-------+--------|
| Fa0 | 5 | 0.05 |
| v0 | 10 | 0.1 |
| V | 66000 | 100 |
| k | 3 | 0.2 |
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.
#+BEGIN_SRC ipython
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
#+END_SRC
#+RESULTS:
:RESULTS:
# Out[67]:
# text/plain
: 0.0050000000000000001
:END:
The uncertainty on Ca is a little trickier. A [[https://en.wikipedia.org/wiki/Propagation_of_uncertainty#Simplification][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 \sigma_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 [[http://kitchingroup.cheme.cmu.edu/blog/2018/10/08/Getting-derivatives-from-implicit-functions-with-autograd/][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.
#+BEGIN_SRC ipython
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}')
#+END_SRC
#+RESULTS:
:RESULTS:
# Out[98]:
# output
: Ca = 0.00500 +\- 0.00016776432898276802
:
:END:
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 [[http://kitchingroup.cheme.cmu.edu/blog/2018/10/11/A-differentiable-ODE-integrator-for-sensitivity-analysis/][differentiable integrator]] will come in handy again!