ODE Sensitivity over time wrt. parameters using cvodes

[Lukas Hebing]

Hello,

After som e reading, I still don`t know, if (and how) it is possible to compute the sensitivity of differential states 'x' to parameters 'p' to over time using either forward- or adjoined sensitivities using cvodes.

The provided example "sensitvity_analysis" describes how to get sensitivities at different time points (e.g. when integrating from t=0 to t=10). 

1) Is there a way to obtain sensitivities at "grid" points or do I have to restart the integrator and simulate from time-point to time-point? Isn`t that computationally more costly? Wouldn`t that lead to different results, depending on the initial state of the sensitivity?

2) Is there a documentation of the commands, which are called in the "sensitvity_analysis" example? For me, it is not clear how the I.factory(...) commands work.

3) Can I "compile" such ODE integrations to a .mex file?

[AY Wer]

Hello,

I want to second the request for an explanation of the I.factory(...) command, if possible.  Especially confusing for me is the string based interface.

Why is it 'adj:qf' in

I_adj = I.factory('I_adj', ['x0', 'z0', 'p', 'adj:qf'], ['adj:x0', 'adj:p'])

and  'adj_qf' in

I_foa = I_adj.factory('I_foa', ['x0', 'z0', 'p', 'adj_qf', 'fwd:p'], ['fwd:adj_x0', 'fwd:adj_p'])

and are the names 'fwd:adj_x0' hardcoded or could we also reference a subset of, e.g., p.

Many thanks.

[AY Wer]

from here I got this example:

from casadi import *


x = MX.sym('x',2); # Two states

p = MX.sym("p")

# Expression for ODE right-hand side
z = p-x[1]**2
rhs = vertcat(z*x[0]-x[1],x[0])

ode = {}         # ODE declaration
ode['x']   = x   # states
ode['p']   = p   # Parameters
ode['ode'] = rhs # right-hand side


# Construct a Function that integrates over 4s
F = integrator('F','cvodes',ode,{'tf':4})

# Start from x=[0;1]
res = F(x0=[0,1],p=1)


print(res["xf"])

# Call with a symbolic parameter
res = F(x0=[0,1],p=p)

# An expression for the integrator sensitivity wrt p
J = jacobian(res["xf"],p)

# A Function to be able to evaluate the expression numerically
JF = Function("JF",[p],[ J ])

print(JF(1))


It would seem that if you are able to use

hessian(res["xf"],p)

but that does not work unfortunately.