How to Update Non-Time-Dependent States in Mixed ODE Problems?

[chooron]

Hello, I’d like to ask for your help. In an ODE problem, after I complete the calculations for a given time step ( t ), I need to update a certain state. This state is not time-dependent but is recalculated based on another variable ( u ), and then used for subsequent calculations.

I have a mixed ODE problem described as follows:

# where s is in m^3 and q_in, q_out are in m^3/s
ds/dt = q_in - q_out

dq_in = new_q_in - q_in 

# where new_q_in is calculated by a function func(q_out), and q_gen is a forcing data.
# We can assume that the function is: func = (q_out, q_gen) -> q_out + q_gen
new_q_in = func(q_out, q_gen)

As you can see, these two balance equations are different: one updates the state ( s ) (in m³), and the other updates the flow rate ( q_{\text{in}} ) (in m³/s, which is time-dependent).

I tried to construct an ODE problem based on this formula. The code is as follows:

using DataInterpolations
using OrdinaryDiffEq

# q_gen
input_arr = [3, 4, 4, 5, 6, 10, 23, 24, 38, 40, 59, 63, 49, 32, 22, 12, 9, 4]
input_itp = LinearInterpolation(input_arr[1,:], collect(1:length(input_arr)))

# ode function
function discharge!(du, u, p, t)
    s_river, q_in = u[1], u[2]
    q_out = (s_river + q_in) / (p[1] + 1)
    du[1] = q_in - q_out
    du[2] = q_out + input_itp(t) - q_in
end

# build problem and solve
prob = DiscreteProblem(discharge!, [0.1, 3], (1, length(input_arr)), [0.3])
sol = solve(prob, FunctionMap())

In the code, I built a DiscreteProblem to solve this problem. However, mathematically, this formula should be continuous, but when constructing an ODEProblem, the results are different:

prob2 = ODEProblem(discharge2!, [0.1, 3], (1, length(input_arr)), [0.3])
sol2 = solve(prob2, Tsit5(), saveat=1.0)

The reason lies in the equation du[2] = q_out + input_itp(t) - q_in, which doesn’t hold in the ODEProblem.

Therefore, I attempted to replace ( u[2] ) by using a DiscreteCallback:

function condition(u, t, integrator)
    return true
end

function affect!(integrator)
    q_in = integrator.u[2]
    q_out = (integrator.uprev[1] + q_in) / (p[1] + 1)
    integrator.u[2] = q_out + input_itp(integrator.t) - q_in
end

cb = DiscreteCallback(condition, affect!)
prob = ODEProblem(discharge!, [0.1, 3], (1, length(input_arr)), [p])
sol = solve(prob, Tsit5(), dt=1.0, callback=cb)
sol_u = Array(sol)

I am not sure if this approach is correct, as the results are still different from earlier (perhaps there’s a mistake in one of the steps).

[janneshb]

Hi,
sorry for my earlier answer, I didn’t get your problem at first and took a closer look now.

Why are you trying to implement this as a DiscreteProblem?
Nothing about your problem seems to be discrete…
In a discrete problem, your function should define a discrete state evolution. I.e. u_{k+1} = f(u_n, p, t_{k+1}).
This is completely different from the ordinary differential equation that you defined, so of course you get very different results from the DiscreteProblem vs the ODEProblem.
Your problem seems to be well defined and I don’t see any issue with taking your first code block and replacing the DiscreteProblem part with

u0 = [0.1, 3]
prob = ODEProblem(discharge!, u0, (1, length(input_arr)), [0.3])
sol = solve(prob, Tsit5())

What you are trying to achieve with the discrete callback? As far as I can tell, you don’t need any discrete dynamics here…

[chooron]

Thank you for your answer.

Indeed, this formula is not a discrete problem, but the issue is that u[2], which is q_in, is not updated like a common state in the ODE.

In the discharge! function,

du[2] = q_out + input_itp(t) - q_in

This calculation formula is theoretically problematic. The unit of q_in is m3/s, so the formula should be:

dq_in/dt = (new_q_in - q_in) /dt

But in the code, dt is ignored, which is why the calculation results differ between DiscreteProblem (dt=1) and ODEProblem(dt!=1).

I also believe that this calculation formula doesn’t need discrete dynamics, but when using DifferentialEquations.jl, I don’t know how to update the state variable q_in. So I considered using a discrete callback, which seems to support updating q_in.

[janneshb]

I honestly don’t really understand your problem. What’s a theoretical sound formulation of this problem, not considering any implementation difficulties?

As far as I can tell, you have an ODE

ds = q_in - q_out

where q_in is computed as

q_in = q_out + q_gen

and q_out is computed as

q_out = (s + q_in) / (p[1] + 1)

q_gen is an input to your system that you compute via linear interpolation of a vector.

Why do you need any discrete dynamics here?

[chooron]

I apologize for not expressing my problem well. The ODE problem you described is somewhat different from my idea:

First, I have an ODE problem:

bash

ds/dt = q_in - q_out # the unit of s is m³, q_in and q_out are m³/s

Then calculate q_out:

q_out = s / (p+1) # here I replace s+q_in with s, because the units of these two variables are different

Finally, I need to obtain the next q_in for calculation based on q_out combined with q_gen, rather than the q_in for this calculation, which is provided by the state u:

q_in(t+dt) = func(q_out(t), q_gen(t))

After obtaining the q_in for the next calculation, I set u[2] = q_in(t+dt) to use it for the next calculation.

I constructed a discrete loop calculation to help you understand this problem:

input_arr = [3 4 4 5 6 10 23 24 38 40 59 63 49 32 22 12 9 4]
s_river = zeros(length(input_arr))
q_in = zeros(length(input_arr))
q_out = zeros(length(input_arr))
s_river0, q_in0 = 0.1, 3
p = [0.3]

for i in eachindex(input_arr)
    q_out[i] = (s_river0) / (p[1] + 1)
    s_river[i] = s_river0 + q_in0 - q_out[i]
    q_in[i] = q_out[i] + input_arr[i]
    s_river0, q_in0 = s_river[i], q_in[i]
    @info [q_out[i] q_in[i] s_river[i]]
end

Then I solved this problem using both ODEProblem and DiscreteProblem, and the results were consistent with the discrete loop calculation, thus confirming your statement:

function discharge2!(du, u, p, t)
    s_river, q_in = u[1], u[2]
    q_out = (s_river) / (p[1] + 1)
    du[1] = q_in - q_out
    du[2] = q_out + input_itp(t) - q_in
end

prob = ODEProblem(discharge2!, [0.1, 3], (1, length(input_arr)), [0.3])
sol3 = solve(prob, Tsit5(), saveat=1.0)

function discharge2(u, p, t)
    s_river, q_in = u[1], u[2]
    q_out = (s_river) / (p[1] + 1)
    [s_river + q_in - q_out, q_out + input_itp(t)]
end

prob = DiscreteProblem(discharge2, [0.1, 3], (1, length(input_arr)), [0.3])
sol4 = solve(prob, FunctionMap(), saveat=1.0)

The results are plotted as follows:

My problem seems to be resolved, but I’m not quite clear why du[2] = q_out + input_itp(t) - q_in can update without considering units. Does du refer to du/dt or just du?

[janneshb]

I don’t really understand the issue you have with units… If your problem is formulated well, all units work out automatically.

But this looks better now:

@chooron

q_in(t+dt) = func(q_out(t), q_gen(t))

However, this

@chooron

du[2] = q_out + input_itp(t) - q_in

doesn’t make a lot of sense, in my opinion. du is a time derivative (in math you would write du/dt), so du[2] would denote the rate of change of q_in, not the next q_in(t+dt). For that you indeed need a discrete problem.

Discrete and continuous problems are completely different. You cannot use the same discharge function for both of them.
Also, u(t + dt) = s_river + q_in - q_out doesn’t seem correct if you’re are assuming sample times other than T = 1.
So, I think the fact that your results make sense are a bit of a coincidence.

For this problem:

@chooron

ds/dt = q_in - q_out

@chooron

q_in(t+dt) = func(q_out(t), q_gen(t))

You do need some hybrid dynamics. I would propose the following:

using DataInterpolations
using OrdinaryDiffEq

input_arr = [3, 4, 4, 5, 6, 10, 23, 24, 38, 40, 59, 63, 49, 32, 22, 12, 9, 4]

function discharge!(du, u, p, t)
    s_river, q_in = u[1], u[2]
    q_out = (s_river + q_in) / (p[1] + 1)
           
    du[1] = q_in - q_out
    du[2] = 0 # q_in is a discrete state, so it will not be changed by the ODE
end

function update_q_in!(integrator)
    q_out = (integrator.u[1] + integrator.u[2]) / (integrator.p[1] + 1)
    integrator.u[2] = q_out + input_arr[Int(round(integrator.t))]
end
cb = PeriodicCallback(update_q_in!, 1.0)
prob = ODEProblem(discharge!, [0.1, 3], (1, length(input_arr)), [0.3])
sol = solve(prob, Tsit5(), dt=1.0, callback=cb)

This will update q_in once a second based on the current q_out and q_gen (from the array).
I don’t understand your problem on a deeper physical level, so I don’t know if this makes sense, but it is a good way to implement such hybrid systems.

Generally, don’t mix continuous and discrete states. If your state evolution is described by an ODE (du[2] = ...) then it probably isn’t also governed by some discrete rule.
Hope this helps

[chooron]

@janneshb

However, this

@chooron

du[2] = q_out + input_itp(t) - q_in

doesn’t make a lot of sense, in my opinion. du is a time derivative (in math you would write du/dt), so du[2] would denote the rate of change of q_in, not the next q_in(t+dt). For that you indeed need a discrete problem.

This is exactly what I wanted to express. du[2] represents the rate of change of q_in, so it cannot be represented by q_out + input_itp(t) - q_in. That’s why I think this might be a discrete problem.

In your code

@janneshb

using DataInterpolations
using OrdinaryDiffEq

input_arr = [3, 4, 4, 5, 6, 10, 23, 24, 38, 40, 59, 63, 49, 32, 22, 12, 9, 4]

function discharge!(du, u, p, t)
    s_river, q_in = u[1], u[2]
    q_out = (s_river + q_in) / (p[1] + 1)
           
    du[1] = q_in - q_out
    du[2] = 0 # q_in is a discrete state, so it will not be changed by the ODE
end

function update_q_in!(integrator)
    q_out = (integrator.u[1] + integrator.u[2]) / (integrator.p[1] + 1)
    integrator.u[2] = q_out + input_arr[Int(round(integrator.t))]
end
cb = PeriodicCallback(update_q_in!, 1.0)
prob = ODEProblem(discharge!, [0.1, 3], (1, length(input_arr)), [0.3])
sol = solve(prob, Tsit5(), dt=1.0, callback=cb)

I’m using DifferentialEquations.jl version 7.14.0, so I’m using DiscreteCallback. The code is as follows:

function discharge!(du, u, p, t)
    s_river, q_in = u[1], u[2]
    q_out = (s_river) / (p[1] + 1)
    du[1] = q_in - q_out
end

function condition(u, t, integrator)
    return true
end

function affect!(integrator)
    q_out = (integrator.u[1]) / (integrator.p[1] + 1)
    integrator.u[2] = q_out + input_itp(integrator.t)
end

cb = DiscreteCallback(condition, affect!)
prob = ODEProblem(discharge!, [0.1, 3], (timeidx[1], timeidx[end]), [0.3])
sol1 = solve(prob, Tsit5(), dt=1.0, callback=cb)

I hope this code should be correct. Thank you for your suggestions！

[janneshb]

@chooron

function condition(u, t, integrator)
    return true
end

This makes no sense, though. If your condition is always true, your discrete time dynamics depend on the step size of the continuous time integrator. That’s why I used a PeriodicCallback with sampling time 1.0. This is crucial. If you solve your ODEProblem with a fixed step size dt=1.0 you are missing out on a lot of nice integration methods provided by DifferentialEquations.jl.

[chooron]

I just find that PeriodicCallback is provided in DiffEqCallbacks.jl, I will try it.

[janneshb]

You can use PeriodicCallback if you use using DifferentialEquations.

The problem with your implementation is, that you did not define a sampling time for your discrete callback. It will always be executed. For that to work you need to define a fixed step size dt = 1.0, like you did.
But this is something you don’t usually want for continuous time systems.
So maybe think about using a periodic callback or adjust your condition() function to only return true once every dt seconds.

[chooron]

I have implemented the method you provided, but it seems to have discrepancies with the calculation results of loop discretization:

# your method
function discharge!(du, u, p, t)
    s_river, q_in = u[1], u[2]
    q_out = (s_river + q_in) / (p[1] + 1)
           
    du[1] = q_in - q_out
    du[2] = 0 # q_in is a discrete state, so it will not be changed by the ODE
end

function update_q_in!(integrator)
    # there may be uprev that we need to use the previous value of s_river
    q_out = (integrator.uprev[1] + integrator.u[2]) / (integrator.p[1] + 1)
    integrator.u[2] = q_out + input_itp(integrator.t)
end
cb = PeriodicCallback(update_q_in!,1.0)
prob = ODEProblem(discharge!, [0.1, 3], (1, length(input_arr)), [0.3])
sol1 = solve(prob, Tsit5(), dt=1.0, callback=cb)

# Use a for loop to calculate
input_arr = [3 4 4 5 6 10 23 24 38 40 59 63 49 32 22 12 9 4]
s_river = zeros(length(input_arr))
q_in = zeros(length(input_arr))
q_out = zeros(length(input_arr))
s_river0, q_in0 = 0.1, 3
p = [0.3]
for i in eachindex(input_arr)
    q_out[i] = (s_river0) / (p[1] + 1)
    s_river[i] = s_river0 + q_in0 - q_out[i]
    q_in[i] = q_out[i] + input_arr[i]
    s_river0, q_in0 = s_river[i], q_in[i]
    @info [q_out[i] q_in[i] s_river[i]]
end

# plot  results
plot(sol1)
plot!(s_river, label="s_river, u[1]")
plot!(q_in, label="q_in, u[2]")

It can be seen that there is a significant discrepancy between these two calculation methods. The problem may lie in the fact that u[2] (q_in) is only updated at dt intervals, so u[2] appears to be in a step-like pattern.

However, in reality, u[2] should be continuously updated based on the constantly changing q_out and q_gen, so the result should be continuous. Perhaps I should consider using other callbacks.

[janneshb]

Ifq_in should be continuously updated then it is in fact not a discrete variable.
I feel like the problem here is not an implementation issue, but a lack of understanding of how the underlying physical system works. Maybe it’s just on my side though.

To summarize: only use du[...] = if you are writing out an ODE, i.e.

\frac{du}{dt} = ...

Use a discrete or periodic callback if you have something like

u_{k+1} = ...

where the you know the time that passes between k and k+1.

If you just have a continuous variable like

q_in = q_out + inout_itp(t)

then you can simply use that in your discharge!function, but q_in should not be a state of the system.

Hope this helps.