Optimizing ODE parameters with constants

New to Julia
,
1

Hi, new to Julia here,
I’m trying to optimize the parameters of an ode system, while keeping some parameter constant.
Using the lotka-volterra as an example,

function f(du,u,p,t)
  du[1] = dx = p[1]*u[1] - u[1]*u[2]
  du[2] = dy = -a*u[2] + u[1]*u[2]
end

u0 = [1.0;1.0]
tspan = (0.0,10.0)
p = [1.5]
prob = ODEProblem(f,u0,tspan,p)

cost_function = build_loss_objective(prob,Tsit5(),L2Loss(t,data),
                                     maxiters=10000,verbose=false)
using Optim
result = optimize(cost_function, 0.0, 10.0)

say I wanted to programmatically pass a constant ‘a’ and then optimize p, then pass a different constant ‘a’ and optimize p again, and so forth for many values of a. I am confused on the syntax for doing this. I was thinking having a function g(a) that includes all the code above and returns the optimized p, but it doesn’t seem to be the best way to do this, especially if I want to try different odes.

Thanks

2

Doing that in a loop would be fine, and throwing @threads over it would parallelize it.

3

Thanks for the quick response (and thanks for the diffeq package in general, as a matlab user it is very refreshing)

If I put both the model f(du,u,p,t) and the optimizer in the same function g(a), then I won’t be able to use f independently of the optimizer. I would instead need to use f (with any given constant a) for other purposes (say plotting for a given p) without running the optimization.
What can I do instead?

4

To qualify a bit what I mean:
I need to be able to do both:

plot f given a, and p

and optimize p given f, and a

5

Closures are your friend here:

function my_dynamics(du,u,p,t,a)

end

f = (du,u,p,t) -> my_dynamics(du,u,p,t,a) # enclose a to now solve with, or enclose other values later
1 Like
6

I rewrote the code based on your suggestion to use closures. I was not familiar with the concept of closures since you mentioned it, so I am not sure that the way I applied is proper, logically and stylistically.
The code seems to work fine, although it’s a little convoluted in my opinion (which may stem from me not using closures properly)

using DifferentialEquations, DiffEqParamEstim
using Optim, Plots

function my_dynamics(du,u,p,t,a)
    du[1] = dx = p[1]*u[1] - u[1]*u[2]
    du[2] = dy = -a*u[2] + u[1]*u[2]
end

f = (du,u,p,t) -> my_dynamics(du,u,p,t,a) # enclose a to now solve with, or enclose other values later

# generate data with a = 1
u0 = [1.0;1.0]
tspan = (0.0,10.0)
p = [1.5]

a = 1
prob = ODEProblem(f,u0,tspan,p)
sol = solve(prob,Tsit5())

t = collect(range(0,stop=10,length=200))
using RecursiveArrayTools # for VectorOfArray
randomized = VectorOfArray([(sol(t[i]) + .01randn(2)) for i in 1:length(t)])
data = convert(Array,randomized)
h = scatter(t,data')

# plot with different values of a
for i = range(1,stop=10)
    new_a = i
    f = (du,u,p,t) -> my_dynamics(du,u,p,t,new_a)
    prob = ODEProblem(f,u0,tspan,p)
    sol = solve(prob,Tsit5())
    plot!(sol)
end

display(h)

# optimize p with given a
function optimization(a)
    f = (du,u,p,t) -> my_dynamics(du,u,p,t,a)
    prob = ODEProblem(f,u0,tspan,p)
    cost_function = build_loss_objective(prob,Tsit5(),L2Loss(t,data),
                   maxiters=10000,verbose=false)

    result = optimize(cost_function, 0.0, 10.0)   
end 

# optimize for a=2
optimization(2)
7

Yeah, in general we’re moving more towards the DiffEqFlux style because it’s a bit more flexible.

https://diffeqflux.sciml.ai/dev/

That then doesn’t have these issues, but it requires you write out a cost function in full.