 # Optimizing ODE parameters with constants

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 = dx = p*u - u*u
du = dy = -a*u + u*u
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

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

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.

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

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

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 = dx = p*u - u*u
du = dy = -a*u + u*u
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)



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.