Why do I get the same parameter after the OptimizationProblem

General Usage
1

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 :slight_smile: 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)

2

What does it print out? What warnings or errors?

3

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]

4

What if you use ADAM?

5

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

6 
using OptimizationOptimisers
OptimizationOptimisers.Adam()
1 Like
7

sorry…
the core code is :slight_smile:
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() .

8

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

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

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:

image
image1402×240 17.9 KB

thanks .

11

Change that to

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

12

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

image
image1117×396 16.2 KB

13

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()

14

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.

16

What is the solve like at the initial parameters?

17

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}:

18

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.

19

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

20

I think @Vaibhavdixit02 was looking at this?

21

@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.