# ODE parameter estimation, parameter not updating

**URL:** <https://discourse.julialang.org/t/ode-parameter-estimation-parameter-not-updating/125783>\
**Category:** General Usage\
**Tags:** optimization, differentialequation\
**Created:** [February 11, 2025, 11:45am UTC](https://discourse.julialang.org/t/ode-parameter-estimation-parameter-not-updating/125783 "2025-02-11T11:45:09Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![Marta\_Pardo\_Araujo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/marta_pardo_araujo/32/208774_2.png) [@Marta\_Pardo\_Araujo](https://discourse.julialang.org/u/Marta_Pardo_Araujo)\
**Post date:** [February 11, 2025, 11:45am UTC](https://discourse.julialang.org/t/ode-parameter-estimation-parameter-not-updating/125783/1 "2025-02-11T11:45:10Z")

</div>

Hi, I have this code to estimate the parameters of a ODE system. I do not know why it is not updating the parameters inside the optimization. I print them and they do not change. This is a test sample of the code:

```julia
import Pkg

using DifferentialEquations, DataFrames, CSV, Plots, LinearAlgebra, ODE, DataInterpolations,
 DiffEqParamEstim, Optimization, Statistics, Dates,ForwardDiff, OptimizationOptimJL, OptimizationBBO, OrdinaryDiffEq,
 OptimizationPolyalgorithms, SciMLSensitivity, Zygote

# Constants
const N = 5
const m_c = 0.001 # probability of mosquito in a car

# Create fake data
const end_ind = 5

# Create fake mobility matrix
const eta = [0.0 7744.15 5240.14 8.57707 5.21319;
        6.9175 0.0 7.23771 45.2212 16.8025;
        61731.9 7.426 0.0 17.7236 54998.8;
        8.40241 44.21 20.1597 0.0 5309.46;
        5830.9 17.5283 5.36841 5226.36 0.0]

# Choose number of patches and IC
pop_init = zeros(end_ind)
pop_init[1] = 1
pop_init[3] = 1

# Non autnomous model ----------------------------------------------------------------------
function fun_na!(du, u, p, t)
  for i in 1:N
      du[i]= sum(j -> p[1] * m_c * eta[i,j] * u[j], 1:N) * (1 - u[i]) -
      p[2]*u[i]
  end
  # println("du:$du[1:5]")
end

# Set initial parameters ------------------------------------------------------------------
t_obs = [592, 957,1323,1688,2053,2418,2784]
t0= t_obs[1] - 10
tf= t_obs[end] + 10
tspan = (t0, tf)
t_vect=1:tf
u0 = pop_init
p = [0.001,0.00000001]

# Create the model
prob = DifferentialEquations.ODEProblem(fun_na!, u0, tspan, p)

# Create the cost function -----------------------------------------------------------------
vec_year = [2007,2008,2009,2010,2011,2012,2013] # Vector of summer times per year

# Function to extract times related to summer per year
function summer_times_by_year(R_M)
    summer_times = Dict{Int, Vector{Int}}()
    for y in vec_year
        summer_times[y] = filter(row -> year(row.date) == y && month(row.date) in 6:8, R_M).time
    end
    return summer_times
end

# Run through RM
summer_t_obs_by_year = Dict(
  2013 => collect(3074:3165),
  2012 => collect(2709:2800),
  2011 => collect(2343:2434),
  2010 => collect(1978:2069),
  2009 => collect(1613:1704),
  2008 => collect(1248:1339),
  2007 => collect(882:973)
)

# Function to compute average probability of occupancy in summer per year and comarca
function average_summer_solution_by_year(sol)
    sol_mat = zeros(Float64,N,length(vec_year))
    j = 1
    for years in vec_year
        times = summer_t_obs_by_year[years]
        vec_sol = Vector{Float64}(undef,N)
        for i in 1:N
            vec_sol[i] = mean(sol(times, idxs = i))
        end  
        sol_mat[:,j] .= vec_sol
        j = j+1 
    end
    return sol_mat
end

# Update p in the ODE
function hanski_prediction(p)
	_prob = remake(prob, p = p)
	sol = solve(_prob, alg_hints=[:stiff]; abstol=1e-8, reltol=1e-8)
  return(sol)
end

# Create matrix observations -------------------------------------------------------------
obs = hanski_prediction(Float64[0.004,0.0000001])
matrix_obs = hcat(obs(t_obs)...)
plot(obs)

# Loss function
function summer_loss_by_year(p)
  println("p:", ForwardDiff.value.(p))
  sol = hanski_prediction(ForwardDiff.value.(p))
  loss = 0.0
  if any((!SciMLBase.successful_retcode(s.retcode) for s in sol)) # Test for not succesfull integration
      print("Loss function INF\n")
      flush(stdout)
      loss = Inf
  else
    summer_avg_by_year = average_summer_solution_by_year(sol)
    loss = sum((summer_avg_by_year .- matrix_obs).^2) 
    println("loss:", ForwardDiff.value.(loss))
    flush(stdout)
  end
  # println("loss:", ForwardDiff.value.(loss))
  # flush(stdout)
  loss
end

# Parameters bounds
low_bound = [0.0,0.0]
up_bound = [2,2]

# Optimizatin -------------------------------------------------------------------------------
adtype = Optimization.AutoForwardDiff() #AutoFiniteDiff()
optf = Optimization.OptimizationFunction((p,x) -> summer_loss_by_year(p), adtype) 

# Loop through different initial parameters parallely --------------------------------------------
init_param = Float64[0.004,0.0000001]
init_param[1] = init_param[1] + (rand(1)/100)[1]
init_param[2] = init_param[2] + (rand(1)/100)[1]

println("init_param:", init_param)

# Create the optimization problem with the new parameters
optprob = Optimization.OptimizationProblem(optf, init_param, lb = low_bound, ub = up_bound) # Create the optimization problem

# Do a try catch to not break when integration problem
try
  optsol = Optimization.solve(optprob,OptimizationOptimJL.LBFGS(), maxiters = 1000)
  print("Finish optimization\n")
catch e
  println("ERROR: $e")
end

```

This is the output:

> p:[0.005178516823542502, 0.0015114340264134889]  
> loss:1.640687614977965  
> p:[0.005178516823542502, 0.0015114340264134889]  
> loss:1.640687614977965  
> p:[0.005178516823542502, 0.0015114340264134889]  
> loss:1.640687614977965  
> retcode: Success  
> u: 2-element Vector{Float64}:  
> 0.005178516823542502  
> 0.0015114340264134889

The parameter p is always the same. But I do not understand why this is happening.

---

<div class="post-metadata">

**Author:** ![ChrisRackauckas](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/chrisrackauckas/32/77_2.png) [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)\
**Post date:** [February 11, 2025, 12:48pm UTC](https://discourse.julialang.org/t/ode-parameter-estimation-parameter-not-updating/125783/3 "2025-02-11T12:48:37Z")

</div>

Don’t manually remove gradient information and then wonder why you get zero gradients. Using the undocumented internals of the AD libraries should only be done with caution because they aren’t supposed to be used by users.

```julia
sol = hanski_prediction(ForwardDiff.value.(p))

```

That’s a very bad idea. If you get rid of derivatives… then the derivative is zero. So the optimization exits.

```julia
function summer_loss_by_year(p)
  println("p:", ForwardDiff.value.(p))
  sol = hanski_prediction(p)
  if any((!SciMLBase.successful_retcode(s.retcode) for s in sol)) # Test for not succesfull integration
      print("Loss function INF\n")
      flush(stdout)
      return Inf
  else
    summer_avg_by_year = average_summer_solution_by_year(sol)
    loss = sum((summer_avg_by_year .- matrix_obs).^2) 
    println("loss:", ForwardDiff.value.(loss))
    flush(stdout)
    return loss
  end
end

```

Then you shouldn’t just assume Float64.

```julia
function average_summer_solution_by_year(sol)    
    stack([[mean(sol(summer_t_obs_by_year[years], idxs = i)) for i in 1:N] for years in vec_year])
end

```

Those two changes then it’s fine.

---

<div class="post-metadata">

**Author:** ![Marta\_Pardo\_Araujo](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/marta_pardo_araujo/32/208774_2.png) [@Marta\_Pardo\_Araujo](https://discourse.julialang.org/u/Marta_Pardo_Araujo)\
**Post date:** [February 11, 2025, 1:16pm UTC](https://discourse.julialang.org/t/ode-parameter-estimation-parameter-not-updating/125783/4 "2025-02-11T13:16:47Z")

</div>

Now it works perfectly. I did not know I was removing gradient information. I did it after the problem with float64 since I thought this was the problem. Thanks a lot for the solution I really appreciate it 🙂