Why do I get the same parameter after the OptimizationProblem

elbert5770 · March 5, 2024, 4:43pm

Here is code that converges. I think there are two problems. One is save_idxs=[1,4] is just going to fail (assuming this means the data is only fit to state variables 1 & 4). State variables 5 & 6 though are sufficient though for convergence. The second problem is that there must be some error in the setup. To troubleshoot, I constructed a loss function without using DiffEqParamEst and it’s fine with Adam, but BFGS with default settings is too aggressive. I also did not use RecursiveArrayTools as it is too easy to mess up the vectorized L2. I just use an array and for loop and it’s so much more transparent.

using DifferentialEquations, Plots, LinearAlgebra, Distributions, Random
using Optimization
# using OptimizationOptimJL
using OptimizationOptimisers
using ForwardDiff

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

function loss(p2,(prob,randomized,num_states,lengtht))
    sol2 =  solve(remake(prob,p=p2), Rodas5P(autodiff = false), adaptive = false, dt = 0.00225, reltol = 1e-12, abstol = 1e-12)
    loss = 0
    
    sol_idxs = [5,6]
    for i in sol_idxs
        sol3 = sol2(t)
        loss += sum((randomized[i,1:lengtht]-(sol3[i,:])).^2)       
    end
    if eltype(p2) == Float64
        @show p2
        @show loss
        # p1 = plot(sol2(t))
        # scatter!(p1,t,[randomized[i,:] for i in 1:7])
        # display(p1) 
    end
    
    return loss
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)
sol = solve(prob, Rodas5P(autodiff = false), adaptive = false, dt = 0.00225, reltol = 1e-12, abstol = 1e-12)
display(plot(sol))
# using RecursiveArrayTools # for VectorOfArray
# using Optimization, ForwardDiff, OptimizationOptimJL
# # using OptimizationOptimisers
t = collect(range(6 * τ, stop = 7 * τ, length = 200))#collect(start,step,stop) #6 * τ, 7 * τ
@show t
randomized = Array{Float64,2}(undef,7,length(t))

for j in eachindex(t)
        randomized[1:7,j] = sol(t[j]) + 0.01*randn(7)
end 


# randomized = VectorOfArray([(sol(t[i]) + 0.01randn(7)) for i in 1:length(t)])

p1 = plot(sol(t))
scatter!(p1,t,[randomized[i,:] for i in 1:7])
display(p1)

p2 = [0.2, 0.3, 0.05, 0.035, 1.31, 1.23, 11.5, 1.6, 0.03]

sol2 =  solve(remake(prob,p=p2), Rodas5P(autodiff = false), adaptive = false, dt = 0.00225, reltol = 1e-12, abstol = 1e-12)
p1 = plot(sol2(t))
scatter!(p1,t,[randomized[i,:] for i in 1:7])
display(p1)


optp = OptimizationFunction(loss,Optimization.AutoForwardDiff())

prob_find_p = OptimizationProblem(optp,p2,(prob,randomized,7,length(t)))
optsol = solve(prob_find_p,OptimizationOptimisers.Adam(0.01), maxiters = 200)
@show optsol.u
sol2 =  solve(remake(prob,p=optsol.u), Rodas5P(autodiff = false), adaptive = false, dt = 0.00225, reltol = 1e-12, abstol = 1e-12)
p1 = plot(sol2(t))
scatter!(p1,t,[randomized[i,:] for i in 1:7])
display(p1)

optsol = solve(remake(prob_find_p,u0=optsol.u),OptimizationOptimisers.Adam(), maxiters = 300)
@show optsol.u
sol2 =  solve(remake(prob,p=optsol.u), Rodas5P(autodiff = false), adaptive = false, dt = 0.00225, reltol = 1e-12, abstol = 1e-12)
p1 = plot(sol2(t))
scatter!(p1,t,[randomized[i,:] for i in 1:7])
display(p1)

optsol = solve(remake(prob_find_p,u0=optsol.u),OptimizationOptimisers.Adam(0.01), maxiters = 500)
@show optsol.u
sol2 =  solve(remake(prob,p=optsol.u), Rodas5P(autodiff = false), adaptive = false, dt = 0.00225, reltol = 1e-12, abstol = 1e-12)
p1 = plot(sol2(t))
scatter!(p1,t,[randomized[i,:] for i in 1:7])
display(p1)

# @show sol(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()

whucsu · March 6, 2024, 4:35am

thanks for your kind heart answers.
firtstly, I used the save_idxs=[1,4], just because the save_idxs=[1] and save_idxs=[4] could be from clinic in the further. maybe I mistaked the use of API,

[quote=“elbert5770, post:22, topic:110243”]
State variables 5 & 6 though are sufficient though for convergence.
[/quote] so what is this mean ?

And I tried your suggestions.

elbert5770:

prob_find_p = OptimizationProblem(optp,p2,(prob,randomized,7,length(t)))
optsol = solve(prob_find_p,OptimizationOptimisers.Adam(0.01), maxiters = 200)
@show optsol.u
sol2 =  solve(remake(prob,p=optsol.u), Rodas5P(autodiff = false), adaptive = false, dt = 0.00225, reltol = 1e-12, abstol = 1e-12)
p1 = plot(sol2(t))
scatter!(p1,t,[randomized[i,:] for i in 1:7])
display(p1)

after this code segment, the Julia gaves exception on my computer like this:

and it seemed that it iterated 6 time, and then crupted.
the code after this has not been test.

elbert5770 · March 6, 2024, 2:52pm

Did you use my entire program or just part? If you ran the entire program, then what Julia version are you using? Operating system? Have you updated packages to the latest version (“up” in package manager, maybe restarting Julia if told to exit, keep running until nothing more is updated). If you copied the whole program, it should converge.

whucsu · March 7, 2024, 1:29am

I update all your codes, after I used big randn weight, it converged, but seem not so stable.
As a julia example, things are solved.
thanks for all.
whu