# Fitting a dynamic system with an exogenous input (nonhomogenous neural ode) via DiffEqFlux
**URL:** https://discourse.julialang.org/t/fitting-a-dynamic-system-with-an-exogenous-input-nonhomogenous-neural-ode-via-diffeqflux/44934
**Category:** New to Julia
**Tags:** package, diffeq
**Created:** [August 14, 2020, 2:19pm UTC](https://discourse.julialang.org/t/fitting-a-dynamic-system-with-an-exogenous-input-nonhomogenous-neural-ode-via-diffeqflux/44934 "2020-08-14T14:19:25Z")
**Posts on this page:** 20
**Page:** 1
### Author:  [@Volker](https://discourse.julialang.org/u/Volker)
#### Post date: [August 14, 2020, 2:19pm UTC](https://discourse.julialang.org/t/fitting-a-dynamic-system-with-an-exogenous-input-nonhomogenous-neural-ode-via-diffeqflux/44934/1 "2020-08-14T14:19:25Z")
Hello,
I´m trying to fit a dynamic system (a hammerstein system, implemented below) with an exogenous input via DiffEqFlux, but actually I don´t know how to put the input excitation into the n-ode-problem. I first tried to follow the example:
[https://diffeqflux.sciml.ai/dev/examples/neural\_ode\_sciml/](https://diffeqflux.sciml.ai/dev/examples/neural_ode_sciml/)
[Example 3: Solving Nonhomogeneous Equations using Parameterized Functions](https://diffeq.sciml.ai/stable/tutorials/ode_example/#Example-3:-Solving-Nonhomogeneous-Equations-using-Parameterized-Functions),
but I wasn´t able to transfer it to my application. I posted below my attempt. Maybe someone could tell me what I´m doing wrong?
```julia
# intialization
using DifferentialEquations, Flux, Optim, DiffEqFlux, DiffEqSensitivity, Plots
# excitation of the system
no_samples = 100
sample_period = 0.1
tsteps = collect(sample_period:sample_period:(no_samples-1)*sample_period)
tspan = (Float32(tsteps[1]), Float32(tsteps[end]))
ex = t->sin(t)
u = ex.(tsteps)
# system
f(x) = (atan(8.0 * x - 4.0) + atan(4.0)) / (2.0 * atan(4.0))
function hammerstein_system(u::Array{Float64})
y= zeros(size(u))
for k in 2:length(u)
y[k] = 0.2 * f(u[k-1]) + 0.8 * y[k-1]
end
return y
end
# system output
y = hammerstein_system(ex.(tsteps))
# model setup
nn_model = FastChain(FastDense(1,50, sigmoid), FastDense(50, 1))
#intial parameter of the model
p_model = initial_params(nn_model)
function predict_neuralode(p)
Array(prob_neuralode(u0, p))
end
function dudt(du, u, p, t)
du[1] = nn_model(ex(t))[1]
end
prob_neuralode = NeuralODE(dudt, tspan, Tsit5(), saveat = tsteps)
prob_neuralode(0.0, p_model)
```
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [August 18, 2020, 3:36pm UTC](https://discourse.julialang.org/t/fitting-a-dynamic-system-with-an-exogenous-input-nonhomogenous-neural-ode-via-diffeqflux/44934/2 "2020-08-18T15:36:11Z")
Optimal control is essentially learning an input excitation, so [https://diffeqflux.sciml.ai/dev/examples/optimal\_control/](https://diffeqflux.sciml.ai/dev/examples/optimal_control/) is an example you might want to look at. Basically, you just pass the function, not much else to it.
---
### Author:  [@Volker](https://discourse.julialang.org/u/Volker)
#### Post date: [August 20, 2020, 1:43pm UTC](https://discourse.julialang.org/t/fitting-a-dynamic-system-with-an-exogenous-input-nonhomogenous-neural-ode-via-diffeqflux/44934/3 "2020-08-20T13:43:14Z")
Thanks. I will try it, but in my point of view the main difference is, that I would like to excitate the system with an excitation signal and in the example the signal is the optimization target. So it is generated by the optimization in contrast to that I would like to specify the signal externally and excitate the system and n-ode with the prepared signal. The reason is that exciatation signals like chirp, aprbs or multisine are useful to gather as much as possible information of the system.
> [@ChrisRackauckas](#):
>
> Basically, you just pass the function, not much else to it.
Okay, I thought, that I would evaluate `ex = t->sin(t)` at the specific time point and so pass the value of it to the nn. I haven’t quite got behind DiffEqFlux’s user interface yet
---
### Author:  [@Volker](https://discourse.julialang.org/u/Volker)
#### Post date: [September 21, 2020, 9:07am UTC](https://discourse.julialang.org/t/fitting-a-dynamic-system-with-an-exogenous-input-nonhomogenous-neural-ode-via-diffeqflux/44934/4 "2020-09-21T09:07:45Z")
I think, I have misunderstood your reply, sorry.
Did you just simply mean this
```julia
prob_neuralode = NeuralODE(dudt, tspan, Tsit5(), saveat = tsteps, ex)
```
with
> [@ChrisRackauckas](#):
>
> Basically, you just pass the function, not much else to it
it still not works, but I think, it evaluates a little bit further. Now I´m getting the following error message.
ERROR: MethodError: Cannot `convert` an object of type var"#17#18" to an object of type Array{Array{Float32,1},1}
Closest candidates are:
Do you maybe know the reason of this error message?
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [September 21, 2020, 12:29pm UTC](https://discourse.julialang.org/t/fitting-a-dynamic-system-with-an-exogenous-input-nonhomogenous-neural-ode-via-diffeqflux/44934/5 "2020-09-21T12:29:47Z")
Could I get what the code looks like? `NeuralODE(dudt, tspan, Tsit5(), saveat = tsteps, ex)` doesn’t make sense because keyword arguments have to go before the keyword arguments.
---
### Author:  [@Volker](https://discourse.julialang.org/u/Volker)
#### Post date: [September 21, 2020, 1:54pm UTC](https://discourse.julialang.org/t/fitting-a-dynamic-system-with-an-exogenous-input-nonhomogenous-neural-ode-via-diffeqflux/44934/6 "2020-09-21T13:54:56Z")
```julia
using DifferentialEquations, Flux, Optim, DiffEqFlux, DiffEqSensitivity, Plots
# excitation of the system
no_samples = 100
datasize = 100
tspan = (0.0f0, 9.9f0)
tsteps = range(tspan[1], tspan[2], length = datasize)
ex = t->sin(t)
u = ex.(tsteps)
f(x) = (atan(8.0 * x - 4.0) + atan(4.0)) / (2.0 * atan(4.0))
function hammerstein_system(u)
y= zeros(size(u))
for k in 2:length(u)
y[k] = 0.2 * f(u[k-1]) + 0.8 * y[k-1]
end
return y
end
y = Float32.(hammerstein_system(ex.(tsteps)))
nn_model = FastChain(FastDense(1,50, sigmoid), FastDense(50, 1))
p_model = initial_params(nn_model)
u0 = Float32.([0.0])
function predict_neuralode(p)
Array(prob_neuralode(u0, p))
end
function dudt(du, u, p, t)
du[1] = nn_model(ex(t))[1]
end
prob_neuralode = NeuralODE(dudt, tspan, Tsit5(), ex, saveat = tsteps)
prob_neuralode(u0, p_model)
```
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [September 21, 2020, 2:39pm UTC](https://discourse.julialang.org/t/fitting-a-dynamic-system-with-an-exogenous-input-nonhomogenous-neural-ode-via-diffeqflux/44934/7 "2020-09-21T14:39:31Z")
The following defines the UDE that you want:
```julia
using DifferentialEquations, Flux, Optim, DiffEqFlux, DiffEqSensitivity, Plots
# excitation of the system
no_samples = 100
datasize = 100
tspan = (0.0f0, 9.9f0)
tsteps = range(tspan[1], tspan[2], length = datasize)
ex = t->sin(t)
u = ex.(tsteps)
f(x) = (atan(8.0 * x - 4.0) + atan(4.0)) / (2.0 * atan(4.0))
function hammerstein_system(u)
y= zeros(size(u))
for k in 2:length(u)
y[k] = 0.2 * f(u[k-1]) + 0.8 * y[k-1]
end
return y
end
y = Float32.(hammerstein_system(ex.(tsteps)))
nn_model = FastChain(FastDense(1,50, sigmoid), FastDense(50, 1))
p_model = initial_params(nn_model)
u0 = Float32.([0.0])
function dudt(du, u, p, t)
du[1] = nn_model(ex(t),p)[1]
end
prob = ODEProblem(dudt,u0,tspan,nothing)
function predict_neuralode(p)
_prob = remake(prob,p=p)
solve(_prob,Tsit5())
end
predict_neuralode(p_model)
```
---
### Author:  [@Volker](https://discourse.julialang.org/u/Volker)
#### Post date: [September 21, 2020, 7:16pm UTC](https://discourse.julialang.org/t/fitting-a-dynamic-system-with-an-exogenous-input-nonhomogenous-neural-ode-via-diffeqflux/44934/8 "2020-09-21T19:16:31Z")
Thank you very much. You helped me a lot
---
### Author:  [@Volker](https://discourse.julialang.org/u/Volker)
#### Post date: [September 22, 2020, 7:25am UTC](https://discourse.julialang.org/t/fitting-a-dynamic-system-with-an-exogenous-input-nonhomogenous-neural-ode-via-diffeqflux/44934/9 "2020-09-22T07:25:46Z")
I have another question. I created the loss function
```julia
#loss(p) = sum(abs2.(y .- [predict_neuralode(p_model).u[i][1] for i in 1:no_samples]))
function loss(p)
sol = predict_neuralode(p)
@show (length(sol))
return sum(abs2.(y[2:end] .- [u[1] for u in sol.u]))
end
res0 = DiffEqFlux.sciml_train(loss,p_model,ADAM(0.005),maxiters=300)
```
and somehow when I evaluate the loss in the main scope `sol` contains 100 elements and inside
```julia
res0 = DiffEqFlux.sciml_train(loss,p_model,ADAM(0.005),maxiters=300)
```
just 99 elements.
Because of this I changed `y` to `y[2:end]` in the loss function, so that I can train, but I don´t know why `sol` has different size.
Additional info:
I increased the datasize to 1000 and the tspan to (0.0, 99.9) and evaluated everything again the problem disappeared, but when I reused 100, it appeared again.
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [September 22, 2020, 8:02am UTC](https://discourse.julialang.org/t/fitting-a-dynamic-system-with-an-exogenous-input-nonhomogenous-neural-ode-via-diffeqflux/44934/10 "2020-09-22T08:02:22Z")
Are you using saveat? What is it like? It might be a single floating point error.
---
### Author:  [@Volker](https://discourse.julialang.org/u/Volker)
#### Post date: [September 22, 2020, 8:11am UTC](https://discourse.julialang.org/t/fitting-a-dynamic-system-with-an-exogenous-input-nonhomogenous-neural-ode-via-diffeqflux/44934/11 "2020-09-22T08:11:58Z")
yes
```julia
function predict_neuralode(p)
_prob = remake(prob,p=p)
solve(_prob,Tsit5(), saveat=0.1)
end
```
and somehow the training isn´t converging. Even the parameters are the same after training.
I checked it like this
Do you have a tip for me?
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [September 22, 2020, 8:19am UTC](https://discourse.julialang.org/t/fitting-a-dynamic-system-with-an-exogenous-input-nonhomogenous-neural-ode-via-diffeqflux/44934/12 "2020-09-22T08:19:47Z")
Make sure the derivative is non-zero and that you’re using the changed `p`. It’s hard to diagnose without code.
As for the saveat issue, I’d like to see exactly how you’re hitting it. My guess is that it’s one of these `0.1+0.1+0.1 != 0.3` issues that needs specific floating point handling somewhere.
---
### Author:  [@Volker](https://discourse.julialang.org/u/Volker)
#### Post date: [September 22, 2020, 9:39am UTC](https://discourse.julialang.org/t/fitting-a-dynamic-system-with-an-exogenous-input-nonhomogenous-neural-ode-via-diffeqflux/44934/13 "2020-09-22T09:39:11Z")
> [@ChrisRackauckas](#):
>
> you’re using the changed `p`
I thought, this is done by the `loss(p)` and using
```julia
DiffEqFlux.sciml_train(loss,p_model,ADAM(0.005),maxiters=300)
```
for the training
> [@ChrisRackauckas](#):
>
> Make sure the derivative is non-zero
You mean the derivative, which I would like to fit?
I think, it is ensured by the excitation of the system with `sin(t)`. If you see the shape of y in the above diagramm, then I would expect \dot{y} would look like similiar then this. (`diff(y, dims=1)`)
Or do you mean something else?
Complete Code:
```julia
using DifferentialEquations, Flux, Optim, DiffEqFlux, DiffEqSensitivity, Plots
plotly()
# excitation of the system
datasize = 100
tspan = (0.0f0, 9.9f0)
tsteps = range(tspan[1], tspan[2], length = datasize)
ex = t->sin(t)
u = ex.(tsteps)
f(x) = (atan(8.0 * x - 4.0) + atan(4.0)) / (2.0 * atan(4.0))
function hammerstein_system(u)
y= zeros(size(u))
for k in 2:length(u)
y[k] = 0.2 * f(u[k-1]) + 0.8 * y[k-1]
end
return y
end
y = Float32.(hammerstein_system(ex.(tsteps)))
nn_model = FastChain(FastDense(1,50, sigmoid), FastDense(50,50, sigmoid), FastDense(50, 1))
p_model = initial_params(nn_model)
u0 = Float32.([0.0])
function dudt(du, u, p, t)
du[1] = nn_model(ex(t),p)[1]
end
prob = ODEProblem(dudt,u0,tspan,nothing)
function predict_neuralode(p)
_prob = remake(prob,p=p)
solve(_prob,Tsit5(), saveat=0.1)
end
function loss(p)
sol = predict_neuralode(p)
return sum(abs2.(y .- [u[1] for u in sol.u]))/length(sol)
end
p_model_ini = copy(p_model)
res0 = DiffEqFlux.sciml_train(loss,p_model,ADAM(),maxiters=1000)
loss(res0.minimizer)
p_model_ini = copy(p_model)
res0 = DiffEqFlux.sciml_train(loss,p_model,ADAM(),maxiters=1000)
res0.minimizer == p_model_ini
res0.minimizer === p_model_ini
y_model = [u[1] for u in predict_neuralode(res0.minimizer).u]
plot([y y_model], label = ["y" "y_model"], ticks=:native,xlabel="time in s", ylabel="y")
```
---
### Author:  [@Volker](https://discourse.julialang.org/u/Volker)
#### Post date: [September 24, 2020, 1:59pm UTC](https://discourse.julialang.org/t/fitting-a-dynamic-system-with-an-exogenous-input-nonhomogenous-neural-ode-via-diffeqflux/44934/14 "2020-09-24T13:59:14Z")
I changed the parameter manually just to see if when I calculate the loss, that the loss will be different. And it is different. So the UDE is reacting on the paramter change. I also looked to your examples to compare the implementations. The example optimal control works fine, but I think the above posted implementation of you is analogously. If I try to debug the code to see whats going inside `DiffEqFlux.sciml_train`, the debugger will stop at line @withprogress progress name=“Training” begin with an error.
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [September 24, 2020, 3:56pm UTC](https://discourse.julialang.org/t/fitting-a-dynamic-system-with-an-exogenous-input-nonhomogenous-neural-ode-via-diffeqflux/44934/15 "2020-09-24T15:56:39Z")
I might get time to take a look later tonight or this weekend.
---
### Author:  [@Volker](https://discourse.julialang.org/u/Volker)
#### Post date: [September 25, 2020, 5:58am UTC](https://discourse.julialang.org/t/fitting-a-dynamic-system-with-an-exogenous-input-nonhomogenous-neural-ode-via-diffeqflux/44934/16 "2020-09-25T05:58:38Z")
Thank you very much for your time.
---
### Author:  [@Volker](https://discourse.julialang.org/u/Volker)
#### Post date: [September 28, 2020, 10:24am UTC](https://discourse.julialang.org/t/fitting-a-dynamic-system-with-an-exogenous-input-nonhomogenous-neural-ode-via-diffeqflux/44934/17 "2020-09-28T10:24:28Z")
I know that I´m maybe a little bit annyoing on this topic, but I´m convinced that UDEs could help me in my task of modeling a stiff system with a data-based approach (this posted problem is just an example and a toy problem to get use to the package). I´m also very pleased for your help and time so far and I can imagine that you have a lot of other things to do. I would like to solve the problem on my own and dig deeper into your implementation, but I found it hard without debugging the code. Is there a general problem with debugging code with macros or is their something wrong with my julia setup,etc. or do I have to do it in a special way (I just use `Juno.@enter DiffEqFlux.sciml_train(loss,p_model,ADAM(0.005), maxiters=100)` in Atom/Juno) ?
Julia Setup:
Julia v1.5.2
DiffEqFlux v1.23.0
DiffEqSensitivity v6.31.6
DifferentialEquations v6.15.0
Flux v0.11.1
Atom v0.12.23
Juno v0.8.4
Debugger v0.6.6
I can evaluated the code, so that the training is successful, but during debugging it stops in the line 117 of train.jl of DiffEqFlux @withprogress progress name=“Training” begin with the error: ERROR: UndefVarError: ProgressLogging not defined. Additionaly I also tried to train via Flux, but also with Flux the parameters aren´t changing, even though they are assigned to be trainable.
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [September 28, 2020, 12:06pm UTC](https://discourse.julialang.org/t/fitting-a-dynamic-system-with-an-exogenous-input-nonhomogenous-neural-ode-via-diffeqflux/44934/18 "2020-09-28T12:06:41Z")
What’s the code you’re running right now?
---
### Author:  [@Volker](https://discourse.julialang.org/u/Volker)
#### Post date: [September 28, 2020, 12:16pm UTC](https://discourse.julialang.org/t/fitting-a-dynamic-system-with-an-exogenous-input-nonhomogenous-neural-ode-via-diffeqflux/44934/19 "2020-09-28T12:16:28Z")
```julia
using DifferentialEquations, Flux, Optim, DiffEqFlux, DiffEqSensitivity, Plots
using LaTeXStrings
plotly()
# excitation of the system
no_samples = 1000
datasize = 1000
tspan = (0.0f0, 99.9f0)
tsteps = range(tspan[1], tspan[2], length = datasize)
ex = t->sin(t)
u = ex.(tsteps)
f(x) = (atan(8.0 * x - 4.0) + atan(4.0)) / (2.0 * atan(4.0))
function hammerstein_system(u)
y= zeros(size(u))
for k in 2:length(u)
y[k] = 0.2 * f(u[k-1]) + 0.8 * y[k-1]
end
return y
end
y = Float32.(hammerstein_system(ex.(tsteps)))
nn_model = FastChain(FastDense(1,50, sigmoid), FastDense(50,50, sigmoid), FastDense(50, 1))
p_model = initial_params(nn_model)
u0 = Float32.([0.0])
function dudt(du, u, p, t)
du[1] = nn_model(ex(t),p)[1]
end
prob = ODEProblem(dudt,u0,tspan,nothing)
# function predict_neuralode(p)
# _prob = remake(prob,p=p)
# solve(_prob,Tsit5(), saveat=0.1)
# end
function predict_neuralode(p)
#_prob = remake(prob,p=p)
solve(prob,Tsit5(), p=p, saveat=0.1)
end
loss(p) = sum(abs2.(y .- Array(predict_neuralode(p))))/length(y)
p_model_ini = copy(p_model)
#p_model_ini = rand(Float32, (size(p_model, 1), 1))
res0 = DiffEqFlux.sciml_train(loss,p_model,ADAM(0.005), maxiters=100)
loss(res0.minimizer)
loss(p_model_ini)
res0.minimizer == p_model_ini
res0.minimizer === p_model_ini
```
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [September 28, 2020, 1:12pm UTC](https://discourse.julialang.org/t/fitting-a-dynamic-system-with-an-exogenous-input-nonhomogenous-neural-ode-via-diffeqflux/44934/20 "2020-09-28T13:12:27Z")
Oh wow… Okay.
```julia
function predict_neuralode(p)
_prob = remake(prob,p=p)
Array(solve(_prob,Tsit5(), saveat=0.1f0))
end
```
That fixes your problem. It’s a bug because
```julia
julia> 99.9 < 99.9f0
true
```
I need to make it robust so that `saveat` as a Number in different precision it will align the endpoint of the adjoint method. Oh man, library writing…
[Next page](https://discourse.julialang.org/t/fitting-a-dynamic-system-with-an-exogenous-input-nonhomogenous-neural-ode-via-diffeqflux/44934.md?page=2)