computationalsciencelaboratory / matlode Goto Github PK

The variable name or Matlode option name is not descriptive or consistent with other names

Hello,
I tried using the adjoint-family integrators to obtain sensitivities of an arbitrary cost functional

w.r.t. interpolation parameters of a linear system

with a linearly interpolated input (using three points) and compared it to the finite difference approximation of the problem. (I attached some commented matlab code as a .zip-file. Use Matlab->Publish->Publish to get a LaTeX-formatted file of the problem and its Jacobians etc.) The sensitivities differ so much that I suspected the finite difference approximation to be wrong. Because of this, I tried implementing a "brute-force" forward sensitivity analysis aswell as an adjoint sensitivity analysis in ode45 by solving the forward sensitivity equations simultaneously and the backward problem after I obtained the forward solution, respectively (this is not included in the minimal working example in the zip file). They both gave the same solution as the finite difference approximation. Interestingly, though, the quadrature value in MATLODE is the same as int_ryp (in provided file) which I calculated numerically using trapz.
In short:

MATLODE output Quadrature after the solver call is Quad = 64.3645 => same as above with trapz!

• Finite difference sensitivity
  

• MATLODE sensitivity
  

=> differ significantly?!

Am I missing something here?
Thank you in advance.

mweMatloteIssue.zip

EDIT: included formated LaTeX here additionally to the matlab comments.

As discussed in previous correspondence with you, passing the Lambda option as a column vector, e.g., to the MATLODE_ERK_ADJ_Integrator function should be the appropriate way. This can also be seen in In "Fatode: A library for forward, adjoint, and tangent linear integration of odes" Appendix C.5, where lambda_N is a column vector, as expected.
This is not working properly, though. The number of adjoint states is initialized using the number of columns returned by the function handle in line 140 with

OPTIONS.NADJ = size(OPTIONS.Lambda(tspan(1), Y0), 2);

So if I provide a column vector, this is always initialized to 1. And, of course, this affects whatever happens in MATLODE_ADJ1_DiscreteIntergrator(...), too.

EDIT: I think I understand now: passing Lambda as a matrix is required if there are multiple functionals to be evaluated, is that correct? Still, I think the dimension check is incorrect?

In the MATLODE_Example_ERK_ADJ_Integrator script, if you change
Time_Interval = [0 20]
to, e.g.,
Time_Interval = [0 10 20]
an error is thrown. I think this is because of
OPTIONS.Quadrature = Quadrature;
in MATLODE_ERK_ADJ_Integrator in line 171 (and possibly others). There, the function handle is overwritten with the numeric integration value after the first for-loop execution. I am not sure how to fix this appropriately.

In the continuous adjoint equations, the gradient of a functional w.r.t. parameters is also depending on how the initial states of the system ODE depend on those parameters. If I understand correctly, this is the option/parameter Y_TLM. But why is this ignored for the adjoint case? Is this somehow implicitly calculated?
Thank you