# Parameter sensitivity of ODE with parameter-dependent event
**URL:** https://discourse.julialang.org/t/parameter-sensitivity-of-ode-with-parameter-dependent-event/103602
**Category:** Numerics
**Tags:** forwarddiff, events, differentialequation
**Created:** [September 7, 2023, 7:25am UTC](https://discourse.julialang.org/t/parameter-sensitivity-of-ode-with-parameter-dependent-event/103602 "2023-09-07T07:25:28Z")
**Posts on this page:** 5
**Page:** 1
### Author:  [@Neodym](https://discourse.julialang.org/u/Neodym)
#### Post date: [September 7, 2023, 7:25am UTC](https://discourse.julialang.org/t/parameter-sensitivity-of-ode-with-parameter-dependent-event/103602/1 "2023-09-07T07:25:28Z")
Hey,
I’m trying to get parameter sensitivities of an ODE problem with a parameter-dependent time-event.
Just throwing AD at the problem doesn’t seem to work:
```julia
using OrdinaryDiffEq
using SciMLSensitivity
using ForwardDiff
# At t = p[2], we assign p[1] <- p[3]
function rhs!(du, u, p, t)
du[1] = -u[1] + p[1]
end
u0 = [1.0]
p_start = [1.2, 2.0, 0.1]
prob = ODEProblem(rhs!, u0, (0.0, 10.0), p_start)
function loss(p)
_prob = remake(prob, p=p)
function condition_disc(u, t, integrator)
return t == integrator.p[2]
end
function condition_cont(u, t, integrator)
return t - integrator.p[2]
end
function affect!(integrator)
# Triggered at t = p[2], use p[3] instead of p[1] for the remaining time
integrator.p[1] = p[3]
end
# sol = solve(_prob, Tsit5(), saveat = 0.0:0.1:10.0, tstops=[ForwardDiff.value(p[2])], callback = DiscreteCallback(condition_disc, affect!), sensealg=ForwardDiffSensitivity())
sol = solve(_prob, Tsit5(), saveat = 0.0:0.1:10.0, tstops=[ForwardDiff.value(p[2])], callback = ContinuousCallback(condition_cont, affect!), sensealg=ForwardDiffSensitivity())
loss = sum(abs2, sol .- 1)
return loss
end
ForwardDiff.gradient(loss, p_start)
```
This yields
```julia
3-element Vector{Float64}:
-2.9178227175008735
0.0
-116.18633113248218
```
But simple finite differences for `p[2]` give
```julia
julia> (loss([1.20, 2.01, 0.1]) - loss([1.2, 1.99, 0.1])) / 0.02
-7.956884892551486
```
Neither `DiscreteCallback` nor `ContinuousCallback` appear to work.
How can I get the derivative wrt. `p[2]` using AD?
---
### Author:  [@akovalov](https://discourse.julialang.org/u/akovalov)
#### Post date: [January 30, 2026, 7:47pm UTC](https://discourse.julialang.org/t/parameter-sensitivity-of-ode-with-parameter-dependent-event/103602/2 "2026-01-30T19:47:25Z")
Hello @Neodym, ~2.5 years after your question I am also looking at the same problem, and finding myself struggling with getting it to work. Have you had any luck on getting the derivatives to materialized with any form of AD? If yes, would you consider sharing the solution.
And generally speaking I was trying to verify the claims/results from [A Comparison of Automatic Differentiation and Continuous Sensitivity Analysis for Derivatives of Differential Equation Solutions](https://arxiv.org/pdf/1812.01892) article about the generalization of the Discrete Local Sensitivity Analysis via AD to hybrid systems, i.e., systems with parametrized events. Unfortunately, I have hard time reproducing the results stated in the Sec. IV for the given system, particularly obtaining the value \frac{\partial y(0)}{\partial a}.
I would be really thankful for any explanation/resource/discussion that could help me verify the results of the paper, and understand the mechanism of injecting the dependence on a into the dataflow of the calculation of y(0).
---
### Author:  [@stevengj](https://discourse.julialang.org/u/stevengj)
#### Post date: [January 30, 2026, 11:38pm UTC](https://discourse.julialang.org/t/parameter-sensitivity-of-ode-with-parameter-dependent-event/103602/3 "2026-01-30T23:38:21Z")
I think this could be one of the many, many problems (once you stray outside of widely used NN components) where AD needs some manual “help”, and it’s well worth learning how to take the derivatives yourself. For ODEs, this is reviewed in [chapter 9 of our course notes](https://arxiv.org/abs/2501.14787) from our [MIT “Matrix Calculus” course](https://github.com/mitmath/matrixcalc/).
Here, you have an ODE of the form
\frac{du}{dt} = f(u,t) = \begin{cases} f\_1(u, t) & t \le t\_0 \\ f\_2(u,t) & t \> t\_0 \end{cases} = f\_1(u,t) + \Theta(t - t\_0) \left(f\_2(u,t) - f\_1(u,t)\right)
where \Theta(t) is the [Heaviside step function](https://en.wikipedia.org/wiki/Heaviside_step_function).
Now, suppose that you want the derivative \frac{\partial u}{\partial t\_0} with respect to the time t\_0 where the jump occurs in the right-hand side. (This is “forward-mode” sensitivity analysis.) This derivative satisfies the linear differential equation:
\frac{d}{dt} \left( \frac{\partial u}{\partial t\_0}\right) = \frac{\partial f}{\partial u} \frac{\partial u}{\partial t\_0} + \frac{\partial f}{\partial t\_0} \\ = \frac{\partial f}{\partial u} \frac{\partial u}{\partial t\_0} - \delta(t - t\_0) \left(f\_2(u(t\_0,t\_0) - f\_1(u(t\_0),t\_0)\right)
Now, you should immediately see why AD will typically have a problem here, unless it is _specifically_ taught how to handle such right-hand-sides: you can probably use AD to compute \partial f/\partial u with no problem, but \partial f/\partial t\_0 yields a [Dirac delta function](https://en.wikipedia.org/wiki/Dirac_delta_function) which AD won’t know what to do with (since it is a distribution rather than an ordinary function).
But if you do it manually, there is no problem with a Dirac delta on the right-hand-side: it just means that \partial u/\partial t\_0 gains a jump discontinuity at t\_0:
\left. \frac{\partial u}{\partial t\_0} \right|\_{t=t\_0^+} = \left. \frac{\partial u}{\partial t\_0} \right|\_{t=t\_0^-} - \left(f\_2(u(t\_0,t\_0) - f\_1(u(t\_0),t\_0)\right)
which you can easily specify via a continuous callback added to the equation \frac{d}{dt} \left( \frac{\partial u}{\partial t\_0}\right) = \frac{\partial f}{\partial u} \frac{\partial u}{\partial t\_0} that holds for t \ne t\_0 (which can be co-evolved with the \frac{du}{dt} ODE).
Note that the initial condition for the sensitivity is \partial u/\partial t\_0 = 0, so this gives a further simplification: the solution of \frac{d}{dt} \left( \frac{\partial u}{\partial t\_0}\right) = \frac{\partial f}{\partial u} \frac{\partial u}{\partial t\_0} is simply \partial u/\partial t\_0 = 0 for t \< t\_0, so you can simply start the solution at t=t\_0 with initial condition \left. \frac{\partial u}{\partial t\_0} \right|\_{t=t\_0^+} = - \left(f\_2(u(t\_0,t\_0) - f\_1(u(t\_0),t\_0)\right).
If you have other parameters besides t\_0, you can differentiate them in the usual way. If you have _lots_ of parameters, and are differentiating a scalar function of the solution (e.g. a loss function), then you may want to implement reverse-mode differentiation. I’ll leave this as an exercise following the description in our course notes, but it is straightforward — again, one term just has a Dirac delta function, which will probably confuse AD but is easy to insert analytically.
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [January 31, 2026, 12:46am UTC](https://discourse.julialang.org/t/parameter-sensitivity-of-ode-with-parameter-dependent-event/103602/4 "2026-01-31T00:46:56Z")
## Reverse Mode
This is already taken into account in the adjoint system. It’s documented here:
> **[Neural Hybrid Differential Equations and Adjoint Sensitivity Analysis | FS](https://frankschae.github.io/post/gsoc-2021/)**
>
> Project summary In this project, we have implemented state-of-the-art sensitivity tools for chaotic dynamical systems, continuous adjoint sensitivity methods for hybrid differential equations, as well as a high level API for automatic...
with a video about it as well:
[ ](https://www.youtube.com/watch?v=iWyWro7T2to)
So if you use SciMLSensitivty.jl it’s all handled. That means any reverse mode that is captured will do it. And there’s tests along these lines for reverse mode here:
>
## Forward Mode
But ForwardDiff.gradient doesn’t capture in the adjoint system. As you can see from the derivation though, the key issue is that you need to differentiate time. This is tested here:
[https://github.com/SciML/OrdinaryDiffEq.jl/blob/v6.106.0/test/ad/autodiff\_events.jl](https://github.com/SciML/OrdinaryDiffEq.jl/blob/v6.106.0/test/ad/autodiff_events.jl)
Now the tricky thing about this is that in order to make direct AD of the solver work in this kind of situation, what needs to happen is that the time span itself must be upgraded to dual numbers, because that’s effectively happening is you need to differentiate through the change in the time point of the event. So it is not sufficient to make `u0` dual valued, you also need to make `tspan`. You can force this by doing something like:
```julia-auto
function loss(p)
_prob = remake(prob, p=p, tspan = eltype(p).(prob.tspan))
```
and that should be the fix that forces differentiation through the callback. Note that this should be happening automatically: the DiffEqBase preprocessing pipeline has a tspan promotion that exists to cover this case:
>
The actual issue here then is a user issue. You can see it in their code:
> [@Neodym](#):
>
> ```julia-auto
> sol = solve(_prob, Tsit5(), saveat = 0.0:0.1:10.0, tstops=[ForwardDiff.value(p[2])], callback = ContinuousCallback(condition_cont, affect!), sensealg=ForwardDiffSensitivity())
>
> ```
`tstops=[ForwardDiff.value(p[2])]` means “drop derivatives on this term”. USERS SHOULD NEVER USE FORWARDDIFF.VALUE OR ONLY DO SO WITH EXTREME CAUTION! It’s not documented for a reason 😅. Now from the analyses above you can see that the differentiation with respect to this tstop value is actually exactly the missing derivative, so setting that dual value to zero is “the bug”. If the user let `tspan` be dual valued and kept the dual here, it would differentiate the callback correctly. So please never use ForwardDiff.value unless you know why you’re doing it.
> [@akovalov](#):
>
> And generally speaking I was trying to verify the claims/results from [A Comparison of Automatic Differentiation and Continuous Sensitivity Analysis for Derivatives of Differential Equation Solutions](https://arxiv.org/pdf/1812.01892) article about the generalization of the Discrete Local Sensitivity Analysis via AD to hybrid systems, i.e., systems with parametrized events. Unfortunately, I have hard time reproducing the results stated in the Sec. IV for the given system, particularly obtaining the value \frac{\partial y(0)}{\partial a}.
The examples there are in the linked package tests. They pass in every release. If you have any trouble reproducing that then open an issue with an MWE.
---
### Author:  [@akovalov](https://discourse.julialang.org/u/akovalov)
#### Post date: [February 2, 2026, 11:28am UTC](https://discourse.julialang.org/t/parameter-sensitivity-of-ode-with-parameter-dependent-event/103602/5 "2026-02-02T11:28:17Z")
Thanks to both of you @ChrisRackauckas, @stevengj for the insight.
Unfortunately I did not manage to find the right test that implements the ODE in the aforementioned manuscript (so if you remember and could kindly point me to it it’d be great!), but I figured out how the event needs to affect the derivative of y in order to incorporate the correct sensitivity.
Just gonna leave it for those here, who might also be interested in a similar question. The system in question is
\begin{align\*} & \dot{x} = f\_x = -a \\ & \dot{y} = f\_y = \begin{cases} & b, x(t) \> 0 \\ & 0 \end{cases}. \end{align\*}
The generic solution for both variables at some time T is as follows:
\begin{align\*} & x(T) = x\_0 + \int\_{0}^{T} f\_x dt = x\_0 + \int\_{0}^{T} (-a) dt \\ & y(T) = y\_0 + \int\_{0}^{T} f\_y dt = y\_0 + \int\_{0}^{t^\*} b dt + \int\_{t^\*}^{T} 0 dt. \end{align\*}
With the resultant derivative wrt a, i.e. \frac{\partial}{\partial a}(\cdot) := \partial\_a (\cdot) (including simplifications)
\begin{align\*} &\partial\_a x(T) = \partial\_a x\_0 - \partial\_a(a \int\_{0}^{T} dt) \\ &\partial\_a y(T) = \partial\_a y\_0 + \underline{\partial\_a \int\_{0}^{t^\*} b dt}, \end{align\*}
where the Leibniz rule of differentiation applies to the underlined term (this is equivalent to dealing with dirac deltas in @stevengj answer), and yields for \partial\_a y(T)
\begin{align\*} &\partial\_a y(T) = b \cdot \partial\_a t^\*, \end{align\*}
where \partial\_a t^\* can be obtained by the invocation of implicit function differentiation upon detection of the event. Since the event is controlled by x\_a(t) = 0, the implicit function theorem tells us that locally t^\*(a) exists which satisfies the x\_a(t) = 0. The implicit derivative wrt a gives
\begin{align\*} & \partial\_a x + \partial\_t x \partial\_a t^\* = 0 \rightarrow \\ & \partial\_a t^\* = -(\partial\_t x)^{-1}\partial\_a x = -(\dot{x})\partial\_a x, \end{align\*}
which gives us the prescription for how to handle the event in forward sensitivity analysis (in pseudo code), i.e.,
```plaintext
y.der = -b * x.der / (-a) = b * x.der / a
```
where x.der term would be picked up by propagating the seeded dual of a.
It is a bit spelled out, but for me personally it was confusing, when I read through the manuscript. Hope somebody else finds it insightful too.