Sensitivities do not match when the running cost explicitly depends on the parameters

[Symplectomorphism]

When the running-cost function, g, depends explicitly on the system parameters, p. I get a mismatch between the gradients reported by adjoint_sensitivities and Calculus.gradient.

I suspect this discrepancy is due to the fact that I calculate dg as the derivative of g w.r.t. only the state, and omit its derivatives w.r.t. the parameters, p.

Below, I provide a minimal example. The dynamics is that of a simple pendulum, controlled at the upward equilibrium point.

An observation: if I make the function g independent of the parameters, then the derivatives do match.

using DifferentialEquations.OrdinaryDiffEq
using DiffEqSensitivity
using DiffEqBase
using Calculus
using ForwardDiff
using QuadGK

function pendulum_eom(dx, x, p, t)
    dx[1] = x[2]
    dx[2] = -sin(x[1]) + (-p[1]*sin(x[1]) + p[2]*x[2])  # Second term is a simple controller that stabilizes π
end

x0 = [0.1, 0.0]
tspan = (0.0, 10.0)
p = [-24.05, -19.137]
​
prob = ODEProblem(pendulum_eom, x0, tspan, p)
sol = solve(prob, Vern9(), abstol=1e-8, reltol=1e-8)

g(x, p, t) = 1.0*(x[1] - π)^2 + 1.0*x[2]^2 + 5.0*(-p[1]*sin(x[1]) + p[2]*x[2])^2
dg(out, y, p, t) = ForwardDiff.gradient!(out, x -> g(x, p, t), y)

res = adjoint_sensitivities(sol,Vern9(),g,nothing,dg,abstol=1e-8,
                                reltol=1e-8,iabstol=1e-8,ireltol=1e-8)

function G(p)
    tmp_prob = remake(prob,p=p)
    sol = solve(tmp_prob,Vern9(),abstol=1e-8,reltol=1e-8)
    res,err = quadgk((t)-> g(sol(t), p, t), 0.0,10.0,atol=1e-8,rtol=1e-8)
    res
end
res2 = Calculus.gradient(G,p)

res' ≈ res2

[ChrisRackauckas]

@YingboMa we need to double check the terms for the continuous adjoint's dg. I think one of the terms might be dg/du instead of dg/dp?

[YingboMa]

The dg/dp term is missing from sensitivity integral.