Why do I get the same parameter after the OptimizationProblem

whucsu · February 15, 2024, 9:02am

Dear All,
I got some questions when using the DiffEqParamEstim.jl

my goal is to get the estimated parameter accroding to the Synthetic Data.
but I got the same parameter as to the firtst “guess” parameter.

the core code is as below:

u0 = [8.0, 8.0, 8.0, 265.0, 0.0, 0.0, 0.0]

p = [0.3, 0.45, 0.06, 0.033, 1.11, 1.13, 11.0, 1.5, 0.03]

tspan = (0, τ * 30)

prob = ODEProblem(Nik_ODE, u0, tspan, p) #ODEProblem(Nik_ODE, u0, tspan, p)

@time sol = solve(prob, Rodas5P(autodiff = false), adaptive = false, dt = 0.00225, reltol = 1e-12, abstol = 1e-12)

using RecursiveArrayTools # for VectorOfArray
using Optimization, ForwardDiff, OptimizationOptimJL, OptimizationBBO,DiffEqParamEstim
t = collect(range(6 * τ, stop = 7 * τ, length = 200))#collect(start,step,stop) #6 * τ, 7 * τ
randomized = VectorOfArray([(sol(t[i]) + 0.01randn(7)) for i in 1:length(t)])
data = convert(Array, randomized)
cost_function = build_loss_objective(prob, Rodas5P(), L2Loss(t, data),
Optimization.AutoForwardDiff(),
maxiters = 100000, verbose = false,save_idxs = [1,4]) #idxs 2 too low ;save_idxs=[1, 2]

#the original p parameter is :[0.3, 0.45, 0.06, 0.033, 1.11, 1.13, 11.0, 1.5, 0.03]

optprob = Optimization.OptimizationProblem(cost_function,[0.1, 0.4, 0.05, 0.023, 1.31, 1.23, 11.5, 1.6, 0.03]) # little change around the p

optsol = solve(optprob,BFGS())#Rodas5P(),BFGS()

the result showed the same as the "[0.1, 0.4, 0.05, 0.023, 1.31, 1.23, 11.5, 1.6, 0.03]

ask for help…
=== the full code is new users can not upload file…sorry
using DifferentialEquations, Plots, LinearAlgebra, Distributions, OffsetArrays, Random, LaTeXStrings

function Valve(R, deltaP)
q = 0.0
if (-deltaP) < 0.0
q = deltaP/R
else
q = 0.0
end
return q

end

function ShiElastance(t, Eₘᵢₙ, Eₘₐₓ, τ, τₑₛ, τₑₚ, Eshift)
τₑₛ = τₑₛτ
τₑₚ = τₑₚτ
#τ = 4/3(τₑₛ+τₑₚ)
tᵢ = rem(t + (1 - Eshift) * τ, τ)

Eₚ = (tᵢ <= τₑₛ) * (1 - cos(tᵢ / τₑₛ * pi)) / 2 +
     (tᵢ > τₑₛ) * (tᵢ <= τₑₚ) * (1 + cos((tᵢ - τₑₛ) / (τₑₚ - τₑₛ) * pi)) / 2 +
     (tᵢ <= τₑₚ) * 0

E = Eₘᵢₙ + (Eₘₐₓ - Eₘᵢₙ) * Eₚ

return E

end

function DShiElastance(t, Eₘᵢₙ, Eₘₐₓ, τ, τₑₛ, τₑₚ, Eshift)

τₑₛ = τₑₛ*τ
τₑₚ = τₑₚ*τ
#τ = 4/3(τₑₛ+τₑₚ)
tᵢ = rem(t + (1 - Eshift) * τ, τ)

DEₚ = (tᵢ <= τₑₛ) * pi / τₑₛ * sin(tᵢ / τₑₛ * pi) / 2 +
      (tᵢ > τₑₛ) * (tᵢ <= τₑₚ) * pi / (τₑₚ - τₑₛ) * sin((τₑₛ - tᵢ) / (τₑₚ - τₑₛ) * pi) / 2
(tᵢ <= τₑₚ) * 0
DE = (Eₘₐₓ - Eₘᵢₙ) * DEₚ

return DE

end

Model parameter values

Eshift = 0.0
Eₘᵢₙ = 0.03
τₑₛ = 0.3
τₑₚ = 0.45
Eₘₐₓ = 1.5
Rmv = 0.06
τ = 1.0

function NIK!(du, u, p, t)
pLV, psa, psv, Vlv, Qav, Qmv, Qs = u
τₑₛ, τₑₚ, Rmv, Zao, Rs, Csa, Csv, Eₘₐₓ, Eₘᵢₙ = p

# 1) Left Ventricle
du[1] = (Qmv - Qav) * ShiElastance(t, Eₘᵢₙ, Eₘₐₓ, τ, τₑₛ, τₑₚ, Eshift) + pLV / ShiElastance(t, Eₘᵢₙ, Eₘₐₓ, τ, τₑₛ, τₑₚ, Eshift) * DShiElastance(t, Eₘᵢₙ, Eₘₐₓ, τ, τₑₛ, τₑₚ, Eshift)
# 2) Systemic arteries 
du[2] = (Qav - Qs ) / Csa     
# 3) Venous
du[3] = (Qs - Qmv) / Csv 
# 4) Left Ventricular Volume
du[4] = Qmv - Qav 
# 5) Aortic Valve flow
du[5]    = Valve(Zao, (pLV - psa)) - Qav
# 6) Mitral Valve flow
du[6]   = Valve(Rmv, (psv - pLV)) - Qmv 
# 7) Systemic flow
du[7]     = (du[2] - du[3]) / Rs
nothing

end

M = [1. 0 0 0 0 0 0
0 1. 0 0 0 0 0
0 0 1. 0 0 0 0
0 0 0 1. 0 0 0
0 0 0 0 0 0 0
0 0 0 0 0 0 0
0 0 0 0 0 0 1. ]

Nik_ODE = ODEFunction(NIK!,mass_matrix=M)
u0 = [8.0, 8.0, 8.0, 265.0, 0.0, 0.0, 0.0]
p = [0.3, 0.45, 0.06, 0.033, 1.11, 1.13, 11.0, 1.5, 0.03]
tspan = (0, τ * 30)
prob = ODEProblem(Nik_ODE, u0, tspan, p) #ODEProblem(Nik_ODE, u0, tspan, p)
@time sol = solve(prob, Rodas5P(autodiff = false), adaptive = false, dt = 0.00225, reltol = 1e-12, abstol = 1e-12)
using RecursiveArrayTools # for VectorOfArray
using Optimization, ForwardDiff, OptimizationOptimJL, OptimizationBBO,DiffEqParamEstim
t = collect(range(6 * τ, stop = 7 * τ, length = 200))#collect(start,step,stop) #6 * τ, 7 * τ
randomized = VectorOfArray([(sol(t[i]) + 0.01randn(7)) for i in 1:length(t)])
data = convert(Array, randomized)
cost_function = build_loss_objective(prob, Rodas5P(), L2Loss(t, data),
Optimization.AutoForwardDiff(),
maxiters = 100000, verbose = false,save_idxs = [1,4]) #idxs 2 too low ;save_idxs=[1, 2]
#the original p parameter is :[0.3, 0.45, 0.06, 0.033, 1.11, 1.13, 11.0, 1.5, 0.03]
optprob = Optimization.OptimizationProblem(cost_function,[0.1, 0.4, 0.05, 0.023, 1.31, 1.23, 11.5, 1.6, 0.03]) # little change around the p
optsol = solve(optprob,BFGS())#Rodas5P(),BFGS()
print(optsol)

ChrisRackauckas · February 15, 2024, 10:13am

What does it print out? What warnings or errors?

whucsu · February 15, 2024, 10:41am

1.238706 seconds (7.12 M allocations: 182.010 MiB, 4.37% gc time, 85.22% compilation time: 98% of which was recompilation)
retcode: Success
u: [0.1, 0.4, 0.05, 0.023, 1.31, 1.23, 11.5, 1.6, 0.03]
Final objective value: Inf

====
the output of u is the same as the initial “guess” p.
of course not as my expection. it should be near to the original p parameter :[0.3, 0.45, 0.06, 0.033, 1.11, 1.13, 11.0, 1.5, 0.03]

ChrisRackauckas · February 15, 2024, 10:43am

What if you use ADAM?

whucsu · February 15, 2024, 10:53am

I am not so familiary with it.
do you mean :
change optprob = Optimization.OptimizationProblem(cost_function,[0.,…
to : optprob = Optimisers.Adam ?
thanks.

ChrisRackauckas · February 15, 2024, 11:09am

using OptimizationOptimisers
OptimizationOptimisers.Adam()

whucsu · February 15, 2024, 11:31am

sorry…
the core code is
cost_function = build_loss_objective(prob, Rodas5P(), L2Loss(t, data),
Optimization.AutoForwardDiff(),
maxiters = 100000, verbose = false,save_idxs = [1,4])
the optimize below got the not expected result
optprob = Optimization.OptimizationProblem(cost_function,[0.1, 0.4, 0.05, 0.023, 1.31, 1.23, 11.5, 1.6, 0.03]) # little change around the p
optsol = solve(optprob,BFGS())#Rodas5P(),BFGS()

I don’t know how to change the above code to Adam() .

whucsu · February 16, 2024, 10:14am

would you please give me the exact suggestion ?
thanks very much.

ChrisRackauckas · February 16, 2024, 11:21am

using OptimizationOptimisers
optsol = solve(optprob,OptimizationOptimisers.Adam())

whucsu · February 17, 2024, 10:37am

the changed code is :
#the original p parameter is :[0.3, 0.45, 0.06, 0.033, 1.11, 1.13, 11.0, 1.5, 0.03]

optprob = Optimization.OptimizationProblem(cost_function,[0.1, 0.4, 0.05, 0.023, 1.31, 1.23, 11.5, 1.6, 0.03]) # little change around the p

optsol = solve(optprob,OptimizationOptimisers.Adam())

I got running error:

thanks .

Vaibhavdixit02 · February 17, 2024, 11:50am

Change that to

optsol = solve(optprob,OptimizationOptimisers.Adam(), maxiters = 100)

whucsu · February 17, 2024, 12:03pm

thanks.
the error was solved.
but after this, I still got the same parameter as the initial “guess” parameter.

the core code is as below:
using RecursiveArrayTools # for VectorOfArray
using Optimization, ForwardDiff, OptimizationOptimJL, OptimizationBBO,DiffEqParamEstim
using OptimizationOptimisers
t = collect(range(6 * τ, stop = 7 * τ, length = 200))#collect(start,step,stop) #6 * τ, 7 * τ
randomized = VectorOfArray([(sol(t[i]) + 0.01randn(7)) for i in 1:length(t)])
data = convert(Array, randomized)
cost_function = build_loss_objective(prob, Rodas5P(), L2Loss(t, data),
Optimization.AutoForwardDiff(),
maxiters = 100000, verbose = false,save_idxs = [1,4])
#the original p parameter is :[0.3, 0.45, 0.06, 0.033, 1.11, 1.13, 11.0, 1.5, 0.03]
optprob = Optimization.OptimizationProblem(cost_function,[0.1, 0.3, 0.05, 0.023, 1.31, 1.23, 11.5, 1.6, 0.03]) # little change around the p
optsol = solve(optprob,OptimizationOptimisers.Adam(),maxiters = 100000)
and the output is still the same as in the optprob

whucsu · February 18, 2024, 7:53am

in the .jl file, after running, the output is like this:
Warning: Instability detected. Aborting
and then gave the same parameter—
Warning: Instability detected. Aborting
└ @ SciMLBase C:\Users\whu.julia\packages\SciMLBase\8XHkk\src\integrator_interface.jl:606
retcode: Success
u: 9-element Vector{Float64}:
0.2
0.3
0.05
0.035
1.31
1.23
11.5
1.6
0.03

maybe the not expected parameter comes from the Instability?

Here I gave the whole code. ask for help…
best wishes.

(new user still can not upload attachment)
using DifferentialEquations, Plots, LinearAlgebra, Distributions, OffsetArrays, Random, LaTeXStrings

function Valve(R, deltaP)
q = 0.0
if (-deltaP) < 0.0
q = deltaP/R
else
q = 0.0
end
return q

end

function ShiElastance(t, Eₘᵢₙ, Eₘₐₓ, τ, τₑₛ, τₑₚ, Eshift)
τₑₛ = τₑₛτ
τₑₚ = τₑₚτ
#τ = 4/3(τₑₛ+τₑₚ)
tᵢ = rem(t + (1 - Eshift) * τ, τ)

Eₚ = (tᵢ <= τₑₛ) * (1 - cos(tᵢ / τₑₛ * pi)) / 2 +
     (tᵢ > τₑₛ) * (tᵢ <= τₑₚ) * (1 + cos((tᵢ - τₑₛ) / (τₑₚ - τₑₛ) * pi)) / 2 +
     (tᵢ <= τₑₚ) * 0

E = Eₘᵢₙ + (Eₘₐₓ - Eₘᵢₙ) * Eₚ

return E

end

function DShiElastance(t, Eₘᵢₙ, Eₘₐₓ, τ, τₑₛ, τₑₚ, Eshift)
τₑₛ = τₑₛτ
τₑₚ = τₑₚτ
#τ = 4/3(τₑₛ+τₑₚ)
tᵢ = rem(t + (1 - Eshift) * τ, τ)

DEₚ = (tᵢ <= τₑₛ) * pi / τₑₛ * sin(tᵢ / τₑₛ * pi) / 2 +
      (tᵢ > τₑₛ) * (tᵢ <= τₑₚ) * pi / (τₑₚ - τₑₛ) * sin((τₑₛ - tᵢ) / (τₑₚ - τₑₛ) * pi) / 2
(tᵢ <= τₑₚ) * 0
DE = (Eₘₐₓ - Eₘᵢₙ) * DEₚ

return DE

end

Model parameter values

Eshift = 0.0
Eₘᵢₙ = 0.03
τₑₛ = 0.3
τₑₚ = 0.45
Eₘₐₓ = 1.5
Rmv = 0.06
τ = 1.0

function NIK!(du, u, p, t)
pLV, psa, psv, Vlv, Qav, Qmv, Qs = u
τₑₛ, τₑₚ, Rmv, Zao, Rs, Csa, Csv, Eₘₐₓ, Eₘᵢₙ = p

# 1) Left Ventricle
du[1] = (Qmv - Qav) * ShiElastance(t, Eₘᵢₙ, Eₘₐₓ, τ, τₑₛ, τₑₚ, Eshift) + pLV / ShiElastance(t, Eₘᵢₙ, Eₘₐₓ, τ, τₑₛ, τₑₚ, Eshift) * DShiElastance(t, Eₘᵢₙ, Eₘₐₓ, τ, τₑₛ, τₑₚ, Eshift)
# 2) Systemic arteries 
du[2] = (Qav - Qs ) / Csa     
# 3) Venous
du[3] = (Qs - Qmv) / Csv 
# 4) Left Ventricular Volume
du[4] = Qmv - Qav 
# 5) Aortic Valve flow
du[5]    = Valve(Zao, (pLV - psa)) - Qav
# 6) Mitral Valve flow
du[6]   = Valve(Rmv, (psv - pLV)) - Qmv 
# 7) Systemic flow
du[7]     = (du[2] - du[3]) / Rs
nothing

end

M = [1. 0 0 0 0 0 0
0 1. 0 0 0 0 0
0 0 1. 0 0 0 0
0 0 0 1. 0 0 0
0 0 0 0 0 0 0
0 0 0 0 0 0 0
0 0 0 0 0 0 1. ]

Nik_ODE = ODEFunction(NIK!,mass_matrix=M)

u0 = [8.0, 8.0, 8.0, 265.0, 0.0, 0.0, 0.0]
p = [0.3, 0.45, 0.06, 0.033, 1.11, 1.13, 11.0, 1.5, 0.03]
tspan = (0, τ * 30)
prob = ODEProblem(Nik_ODE, u0, tspan, p) #ODEProblem(Nik_ODE, u0, tspan, p)
@time sol = solve(prob, Rodas5P(autodiff = false), adaptive = false, dt = 0.00225, reltol = 1e-12, abstol = 1e-12)

using RecursiveArrayTools # for VectorOfArray
using Optimization, ForwardDiff, OptimizationOptimJL, OptimizationBBO,DiffEqParamEstim
using OptimizationOptimisers
t = collect(range(6 * τ, stop = 7 * τ, length = 200))#collect(start,step,stop) #6 * τ, 7 * τ
randomized = VectorOfArray([(sol(t[i]) + 0.01randn(7)) for i in 1:length(t)])
data = convert(Array, randomized)
cost_function = build_loss_objective(prob, Rodas5P(), L2Loss(t, data),
Optimization.AutoForwardDiff(),tspan=(6 * τ, 7 * τ),save_idxs = [1,4]) #idxs 2 too low ;save_idxs=[1, 2]
#the original p parameter is :[0.3, 0.45, 0.06, 0.033, 1.11, 1.13, 11.0, 1.5, 0.03]
optprob = Optimization.OptimizationProblem(cost_function,[0.2, 0.3, 0.05, 0.035, 1.31, 1.23, 11.5, 1.6, 0.03]) # little change around the p
optsol = solve(optprob,Optim.BFGS())
#optsol = solve(optprob,OptimizationOptimisers.Adam(),maxiters = 100000)#OptimizationOptimisers.Adam(),BFGS()

ChrisRackauckas · February 18, 2024, 1:49pm

Yes, if your differential equation is unstable at the chosen starting parameters and requires aborting then it will not be able to update the parameters. That’s the reason for the warning.

ChrisRackauckas · February 22, 2024, 4:26am

What is the solve like at the initial parameters?

whucsu · February 23, 2024, 5:45am

it seemed that the initial parameters are good.
the output is like this:

retcode: Success
Interpolation: specialized 4rd order “free” stiffness-aware interpolation
t: 13335-element Vector{Float64}:

whucsu · February 24, 2024, 6:27am

I checked the initial parameters of p=[0.25, 0.4, 0.05, 0.023, 1.31, 1.23, 11.5, 1.6, 0.03] again, the return code is retcode: Success.
is there any method that could print out the result when the system is instability while updating the parameters ?
but it just seemed that the parameters have never been updated .

thanks.

whucsu · February 28, 2024, 3:54am

is there any other method that could be used to check the stability ?
debug ? check the level of output ?
thanks. best wishes.
whu

ChrisRackauckas · March 5, 2024, 9:02am

I think @Vaibhavdixit02 was looking at this?

Vaibhavdixit02 · March 5, 2024, 4:32pm

@whucsu there are two things off in your script that make it not work. The way to diagnose it was run the cost_function(p) which gives Inf in your current script.
I am pasting the lines that were changed to make it work

data = convert(Array, randomized)[[1,4], :]
cost_function = build_loss_objective(prob, Rodas5P(autodiff = false), L2Loss(t, data),
                                    Optimization.AutoForwardDiff(),save_idxs = [1,4])

Basically, since you want to use save_idxs you should only pass those states’ data. Also the Rodas5P didn’t have the autodiff = false in the cost_function definition.

Topic		Replies	Views
Discrepancy in DiffEqParamEstim Optimization (Mathematical)	1	366	December 24, 2022
Issue with example in DiffEqParamEstim.jl New to Julia diffeq , examples , optim , differentialequation	2	458	June 21, 2022
Optimization on Vector-ODEs for parameter estimation New to Julia question	3	382	August 21, 2023
DiffEqParamEstim.jl and multiple_shooting_objective General Usage diffeq , optim	6	1134	May 17, 2021
Optimizing a 2 parameter simple ODE Optimization (Mathematical)	1	308	August 25, 2022

Why do I get the same parameter after the OptimizationProblem

Model parameter values

Model parameter values

Related topics