# ModelingToolkit.jl ODESystems with a forcing function: how to autodiff?

**URL:** <https://discourse.julialang.org/t/modelingtoolkit-jl-odesystems-with-a-forcing-function-how-to-autodiff/58823>\
**Category:** Modelling & Simulations\
**Created:** [April 8, 2021, 12:29pm UTC](https://discourse.julialang.org/t/modelingtoolkit-jl-odesystems-with-a-forcing-function-how-to-autodiff/58823 "2021-04-08T12:29:41Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![Dhruva2](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dhruva2/32/18475_2.png) [@Dhruva2](https://discourse.julialang.org/u/Dhruva2)\
**Post date:** [April 8, 2021, 12:29pm UTC](https://discourse.julialang.org/t/modelingtoolkit-jl-odesystems-with-a-forcing-function-how-to-autodiff/58823/1 "2021-04-08T12:29:41Z")

</div>

Hi!

Lots of differential equations models incorporate an external input function (e.g. a control system doing reference tracking). When I build models of this form in ModelingToolkit.jl, they don’t seem to be differentiable. What am I doing wrong?

MWE:

A linear-time-invariant system being forced by a sine wave input. the A matrix constitutes the parameters of the ODE.

```julia
using ModelingToolkit, Zygote, OrdinaryDiffEq, DiffEqSensitivity, LinearAlgebra

function simple_sys() #u is input
    @variables t x[1:2](t) u(t)
    @parameters p[1:4]
    D = Differential(t)
    A = reshape(p,2,2)
    B = reshape([2.,3], 2, 1)

    eqs = D.(x) .~ A*x .+ B*u
    defaults = Dict(vcat(p .=> [1., -2., 3., -0.5], x .=> ones(2), u => 0.))
    return ODESystem(eqs; name=:lp, defaults)
end
od = simple_sys()

function connect_to_input(od, input; name)
    @variables t f(t)
    eqs = [
        f ~ input(t),
        f ~ od.u
    ]
    ODESystem(eqs; systems = [od], name = name, defaults = Dict(f .=> input(0.)))
end

@named con = connect_to_input(od, sin)
final = structural_simplify(con)
prob = ODEProblem(final, [], (0., 10.),[])
tsteps = 0:0.5:10
nom_sol = solve(prob, Tsit5(), saveat = tsteps, sensealg = InterpolatingAdjoint())

function loss(p)
    nprob = remake(prob; p=p)
    sol = solve(nprob, Tsit5(), saveat = tsteps, sensealg = InterpolatingAdjoint())
    return sum(abs2, sol[Symbol("lp₊x₁(t)")] .- sin.(tsteps))
end

Zygote.gradient(loss, prob.p)

```

Stacktrace is:

```julia
ERROR: Can't differentiate foreigncall expression
Stacktrace:
 [1] error(s::String)
   @ Base ./error.jl:33
 [2] Pullback
   @ ./boot.jl:480 [inlined]
 [3] (::typeof(∂(Symbol)))(Δ::Nothing)
   @ Zygote ~/.julia/packages/Zygote/lwmfx/src/compiler/interface2.jl:0
 [4] Pullback
   @ ~/.julia/dev/FeedbackLearning/scripts/undifferentiable.jl:35 [inlined]
 [5] (::typeof(∂(loss)))(Δ::Float64)
   @ Zygote ~/.julia/packages/Zygote/lwmfx/src/compiler/interface2.jl:0
 [6] (::Zygote.var"#41#42"{typeof(∂(loss))})(Δ::Float64)
   @ Zygote ~/.julia/packages/Zygote/lwmfx/src/compiler/interface.jl:41
 [7] gradient(f::Function, args::Vector{Float32})
   @ Zygote ~/.julia/packages/Zygote/lwmfx/src/compiler/interface.jl:59
 [8] top-level scope
   @ REPL[97]:1

```

How should I set up a forcing function so that autodiff works?

Thanks a lot!

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [April 22, 2021, 7:54am UTC](https://discourse.julialang.org/t/modelingtoolkit-jl-odesystems-with-a-forcing-function-how-to-autodiff/58823/2 "2021-04-22T07:54:31Z")

</div>

The problem was:

```julia
function loss(p)
    nprob = remake(prob; p=p)
    sol = solve(nprob, Tsit5(), saveat = tsteps, sensealg = InterpolatingAdjoint())
    return sum(abs2, sol[Symbol("lp₊x₁(t)")] .- sin.(tsteps))
end

Zygote.gradient(loss, prob.p)

# Works
x = Symbol("lp₊x₁(t)")
function loss(p)
    nprob = remake(prob; p=p)
    sol = solve(nprob, Tsit5(), saveat = tsteps, sensealg = InterpolatingAdjoint())
    return sum(abs2, sol[x] .- sin.(tsteps))
end

Zygote.gradient(loss, prob.p)

```

which means it’s fixed by `Zygote.@nograd Symbol`. Thus [Add nograd to Symbol constructor by ChrisRackauckas · Pull Request #952 · FluxML/Zygote.jl · GitHub](https://github.com/FluxML/Zygote.jl/pull/952) is the true fix.

---

<div class="post-metadata">

**Author:** ![Dhruva2](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dhruva2/32/18475_2.png) [@Dhruva2](https://discourse.julialang.org/u/Dhruva2)\
**Post date:** [April 28, 2021, 3:54pm UTC](https://discourse.julialang.org/t/modelingtoolkit-jl-odesystems-with-a-forcing-function-how-to-autodiff/58823/3 "2021-04-28T15:54:08Z")

</div>

Thanks, yes that was indeed the problem, which is now solved.