* Sensitivity analysis with odeint and autograd
:PROPERTIES:
:categories: autograd,ode
:date: 2019/09/13 09:56:09
:updated: 2019/09/13 09:56:09
:org-url: http://kitchingroup.cheme.cmu.edu/org/2019/09/13/Sensitivity-analysis-with-odeint-and-autograd.org
:permalink: http://kitchingroup.cheme.cmu.edu/blog/2019/09/13/Sensitivity-analysis-with-odeint-and-autograd/index.html
:END:
In this [[http://kitchingroup.cheme.cmu.edu/blog/2018/10/11/A-differentiable-ODE-integrator-for-sensitivity-analysis/][previous post]] I showed a way to do sensitivity analysis of the solution of a differential equation to parameters in the equation using autograd. The basic approach was to write a differentiable integrator, and then use it in a function so that autograd could take the derivative.
Since that time, autograd has added [[https://github.com/HIPS/autograd/blob/master/autograd/scipy/integrate.py][derivative support]] for =scipy.integrate.odeint=. In this post we examine that. As usual with autograd, we have to import the autograd version of numpy, and the autograd version of odeint. We will find the derivative of the solution to an ODE (which is an array) so we need to also import the jacobian function. Finally, there is a subtle, and non-obvious requirement that we need to import the autograd tuple. That ensures that the variables are differentiable through the tuple we will use for the arguments.
The differential equation we solve returns the concentration of a species as a function of time, and the solution depends on two parameters, i.e. $C = f(t; k_1, k_{-1})$, and we are interested in the time-dependent sensitivity of $C$ with respect to those parameters. The approach we use is to define a function that has those parameters as arguments. The function will solve the ODE and return the time-dependent solution. First we make that solution, mostly to see that the autograd version of odeint works.
#+BEGIN_SRC ipython
import autograd.numpy as np
from autograd.scipy.integrate import odeint
from autograd import jacobian
from autograd.builtins import tuple
import matplotlib.pyplot as plt
Ca0 = 1.0
k1 = k_1 = 3.0
tspan = np.linspace(0, 0.5)
def C(K):
k1, k_1 = K
def dCdt(Ca, t, k1, k_1):
return -k1 * Ca + k_1 * (Ca0 - Ca)
sol = odeint(dCdt, Ca0, tspan, tuple((k1, k_1)))
return sol
plt.plot(tspan, C([k1, k_1]))
plt.xlim([tspan.min(), tspan.max()])
plt.xlabel('t')
plt.ylabel('C');
#+END_SRC
#+RESULTS:
:results:
# Out [50]:
# text/plain
: