casadi / casadi

CasADi is a symbolic framework for numeric optimization implementing automatic differentiation in forward and reverse modes on sparse matrix-valued computational graphs. It supports self-contained C-code generation and interfaces state-of-the-art codes such as SUNDIALS, IPOPT etc. It can be used from C++, Python or Matlab/Octave.
http://casadi.org
GNU Lesser General Public License v3.0
1.67k stars 370 forks source link

Inconsistent behavior of ode sensitivity analysis in CasADi #3783

Open gaowutong opened 1 month ago

gaowutong commented 1 month ago

The code demonstrates an issue with the calculation of the Jacobian of a subset of the output vector xf. Specifically, the calculation works fine for the full sensitivity matrix and the first 3 rows, but fails for the first 2 rows.

from casadi import *

r = SX.sym('r', 3)
v = SX.sym('v', 3)
g = -r / norm_2(r) ** 3
x = vertcat(r, v)
xdot = vertcat(v, g)

dae = {'x': x, 'ode': xdot}
I = integrator('I', 'cvodes', dae, 0, 1)

x0 = MX.sym('x0', 6)
xf = I(x0=x0)['xf']

x0_eval = [1, 0, 0, 0, 1, 0]
print(Function('dxf_dx0', [x0], [jacobian(xf, x0)])(x0_eval))  # full sensitivity matrix, works fine
print(Function('dxf_dx0', [x0], [jacobian(xf[:3], x0)])(x0_eval)) # first 3 rows, works fine
print(Function('dxf_dx0', [x0], [jacobian(xf[:2], x0)])(x0_eval)) # first 2 rows, fails

The first and second output is:

[[1.89694, 0.386823, 4.79716e-24, 1.22829, 0.188871, -4.59979e-23], 
 [0.516683, 0.751624, 4.6718e-24, 0.211322, 0.971331, -5.74475e-24], 
 [1.19833e-22, 3.37775e-22, 0.540302, 2.566e-22, 3.18579e-22, 0.841471], 
 [1.55308, 0.956449, 2.08845e-23, 1.49675, 0.643786, -6.14408e-24], 
 [1.56796, -0.0678248, 1.57219e-23, 0.773646, 1.15182, -2.03127e-23], 
 [-1.64163e-23, -8.78783e-23, -0.841471, 6.24318e-23, -1.2098e-22, 0.540302]]

[[1.89694, 0.386823, 4.79716e-24, 1.22829, 0.188871, -4.59979e-23], 
 [0.516683, 0.751624, 4.6718e-24, 0.211322, 0.971331, -5.74475e-24], 
 [1.19833e-22, 3.37775e-22, 0.540302, 2.566e-22, 3.18579e-22, 0.841471]]

While the third got unexpected error:

rhsQB failed: .../casadi/core/oracle_function.cpp:514: Assertion "it!=all_functions_.end()" failed:
No function "quadB" in asens2_I. Available functions: daeB,daeF,jacF.
At t = 1, the quadrature right-hand side routine failed in an unrecoverable manner.
Error occured while integrating backward problem # 0
Function asens2_I (0x55630bcd9190)
Input 0 (x0): [1, 0, 0, 0, 1, 0]
Input 1 (z0): 0x1
Input 2 (p): 0x1
Input 3 (u): 0x1
Input 4 (adj_xf): 
[[1, 0], 
 [0, 1], 
 [0, 0], 
 [0, 0], 
 [0, 0], 
 [0, 0]]
Input 5 (adj_zf): 0x2
Input 6 (adj_qf): 0x2
Function adj2_I (0x55630be01cd0)
Input 0 (x0): [1, 0, 0, 0, 1, 0]
Input 1 (z0): 0x1
Input 2 (p): 0x1
Input 3 (u): 0x1
Input 4 (adj_xf): 0x0
Input 5 (adj_zf): 0x0
...
Input 17 (adj2_adj_x0): 0x0
Input 18 (adj2_adj_z0): 0x0
Input 19 (adj2_adj_p): 0x0
Input 20 (adj2_adj_u): 0x0

By the way, the code works fine for 'rk' and 'collocation' integrator, but fails for 'cvodes' and 'idas'.

tmmsartor commented 1 month ago

Just to add a bit on this, the problem seems that between 3 to 2 rows it switches from using Forward to Adjoint. In case of cvodes and idas some required functions (quadB) is not created when parameter size is 0, in fact using dae = {'x': x, 'p': SX.sym('p',1), 'ode': xdot} does not expose the issue.

gaowutong commented 1 week ago

Thanks, this is a way to solve the issue.

tmmsartor commented 1 week ago

I am happy the workaround is enough for you but I think this issue is not completely solved, I keep it open for now.