Trouble with Optimization.jl errors

Specific Domains Optimization (Mathematical)
, ,
1

Hello all. I’m running into an issue with Optimization.jl when trying to estimate parameters of ODEs. The specific issue is that some optimizers work while others don’t. Furthermore, different optimizers throw different errors. This is not an issue of the optimizers failing to converge, but rather they throw errors before starting.

The following is a simple example. Here a single ODE is setup with some synthetic data. Much of this follows online SciML tutorials. When trying to fit the ODE to the data, PolyOpt works but CMAES and SimulatedAnnealing do not. I have not tried an exhaustive search to see what works and what doesn’t. Note that I’m using ComponentArrays for my parameter vectors, but this is not the cause. If I use standard vectors I get similar failurs (though different errors in some cases).

Does anyone have any guidance for approaching this issue? I work with systems of 10-50 ODEs and need the ability to use both gradient based (PolyOpt is perfect) and particle methods (Genetic Algo, etc.).

I am trying to determine the best way forward here. I don’t know the internals of this ecosystem well so any guidance would be helpful. 1) Am I doing something silly in the code below that is causing problems with my interface to the optimizers? 2) Is there just a better way to use julia’s optimization infrastructure with ODE optimization?


using Pkg
Pkg.activate("SciML_ParamEst_Env/") # Ignore of course.

using OrdinaryDiffEq,
      Optimization, OptimizationPolyalgorithms, OptimizationEvolutionary,
      OptimizationOptimJL

using ForwardDiff
using ComponentArrays

##### Define the logistic growth model #####
function logistic_rhs!(du, u, p, t)
    x = u[1]
    r, k = p.Params
    du[1] =  r * x * (1.0 - x/k)
    
    return nothing
end

# Generating Parameters
u0, r, k = 0.2, 1.3, 2.0;
p_true = ComponentArray( IC = [u0], Params = [r, k] )

# Simulation interval and intermediary points
tspan = (0.0, 5.0);
tsteps = 0.0:0.1:5.0;

# Setup the ODE problem and simulate some synthetic data
base_prob = ODEProblem(logistic_rhs!, p_true.IC, tspan, p_true);
sol_exact = solve(base_prob, Tsit5(),saveat = tsteps);

# Define SSE loss function.
# Note that 'p' pulls in some things we want to avoid passing as globals
function loss(new_params , p) 
    
    odeprob = p[1] # ODEProblem stored as parameters to avoid using global variables
    tsteps = p[2]; # time steps to compare the solution
    data = p[3]; # data to compare the solution with
    
    sol = solve(odeprob, Tsit5() , u0 = new_params.IC , 
                                    p = new_params, 
                                    saveat = tsteps)
    data_sim = Array(sol)

    return sum((data .- data_sim) .^ 2) / length(data)
end

##### Perform optimization #####

# This stores the basic odeprob structure to avoid repeaded generation.
# Also stores the exact solution (and time points) to avoid globals.
params_fixed = [base_prob, tsteps, Array(sol_exact)];

# Initial guess
param_start = ComponentArray(IC = [1.0], Params = [1.0, 1.0] );

adtype = Optimization.AutoForwardDiff()
optf = Optimization.OptimizationFunction((x, p) -> loss(x,params_fixed), adtype)
optprob = Optimization.OptimizationProblem(optf, param_start)

# WORKS
result_ode = Optimization.solve(optprob, PolyOpt(),
    maxiters = 200)

# ERROR
# "MethodError: no method matching value!"
result_ode = Optimization.solve(optprob, CMAES(μ = 40, λ = 100),
    maxiters = 200)

# ERROR
# DimensionMismatch: arrays could not be broadcast to a common size
result_ode = Optimization.solve(optprob, SimulatedAnnealing(),
    maxiters = 200)
2

I don’t know much about Optimization, but these just look like bugs. Can you post the full stack trace of the errors?

You might want to post on GitHub · Where software is built, otherwise @Vaibhavdixit02 is the one to ask.

1 Like