Open elbert5770 opened 1 year ago
@avik-pal could this be due to the new termination conditions?
That is not entirely surprising to me. The tolerances for time-stepping specified via solve
but tolerances for termination is specified via abstol
and reltol
in DynamicSS
https://github.com/SciML/SteadyStateDiffEq.jl/blob/063e95843122abc5a983845d9b2cd887e6981f75/src/algorithms.jl#L47
The abstol and reltol are inside DynamicSS, not attached to solve. Moving them outside of DynamicSS into solve makes the termination much worse, as expected. But these are definitely inside DynamicSS and affecting termination.
Corrected in edit: xxx[Okay, I think the root issues are the aggressive time step growth and the reltol = 1e-6 is too large. The delta t is growing fast and the residuals are bouncing around. Setting dtmax in solve and lowering reltol to say 1e-10 leads to a well-behaved algorithm for both the FiniteDiff and ForwardDiff gradient. sol = solve(remake(probss,u0=eltype(kf_guess).(u0),p=(p_guess,p_fix,kf_guess,plateau,kd,start_time,end_time)),DynamicSS(AutoTsit5(Rosenbrock23()),reltol=1e-10),dtmax=1.0)
] xxx
Nope, this only worked because I had changed p_guess to [10.0,3.0]. With the original p_guess, it terminates too quickly with reltol=1e-10.
As Avik mentioned, the time stepping can be affected by the abstol and reltol to solve, so this might be better than dtmax:
sol = solve(remake(probss,u0=eltype(kf_guess).(u0),p=(p_guess,p_fix,kf_guess,plateau,kd,start_time,end_time)),DynamicSS(AutoTsit5(Rosenbrock23()),abstol=1e-16,reltol=1e-16,tspan=1e6),abstol=1e-11,reltol=1e-11)
Here, unfortunately, the solver must run to t=1e6 to get good results and avoid early termination consistently, or hit these very low values of reltol and abstol within DynamicSS.
Is this fixed with the new termination conditions? https://docs.sciml.ai/NonlinearSolve/stable/basics/TerminationCondition/
It appears to be exactly the same, but I'm only 99% sure of my implementation of TerminationCondition:
using FiniteDiff
using ForwardDiff
using Optimization
using OptimizationOptimJL
using NonlinearSolve
using DifferentialEquations
using SteadyStateDiffEq
using Plots
function sys!(du, u, (p,p_fix,kf,plateau,kd,start_time,end_time), t)
#@show t,forcing_function(t,plateau,kd,start_time,end_time),p,p_fix,kf
du[1] = kf[1]*(1.0 - forcing_function(t,plateau,kd,start_time,end_time)) - p[1]*u[1]
du[2] = p[1]*u[1] - p_fix[1]*u[2]
du[3] = p_fix[1]*u[2] - p[2]*u[3]
if t > 50
# @show u
# @show du
@show t
end
return nothing
end
function forcing_function(t,plateau,kd,start_time,end_time)
if t < start_time
return plateau*t
else
if t <= end_time
return plateau
else
if t > end_time && t < Inf
return plateau*(exp(-kd*(t-end_time)))
else
return 0.0
end
end
end
end
function set_kf_outer(probss,u03,termination_condition)
function set_kf_inner(kf_guess,(u0,p_guess,p_fix,plateau,kd,start_time,end_time))
# Solve for steady state with current p_guess and an initial guess of kf and u0
plateau = 0.0
#sol = solve(remake(probss,u0=eltype(kf_guess).(u0),p=(p_guess,p_fix,kf_guess,plateau,kd,start_time,end_time)),DynamicSS(AutoTsit5(Rosenbrock23()),abstol = 1e-14,reltol = 1e-14))
sol = solve(remake(probss,u0=eltype(kf_guess).(u0),p=(p_guess,p_fix,kf_guess,plateau,kd,start_time,end_time)),DynamicSS(AutoTsit5(Rosenbrock23()),termination_condition=termination_condition))
# Compare kf at steady state to measured u0[3]
# @show sol.u
# @show sol.resid
# @show sol.prob
@show (sol[3]-u03)^2
return (sol[3]-u03)^2
end
end
function main()
#### Parameters and state variables
kf0= [10.0] #true value of zero order production rate
p0 = [1.0,0.25,0.5] #true values of parameters
p_fix = [p0[2]] # a parameter that is known, note that the system is otherwise unidentifiable
p_var = [p0[1],p0[3]] # Unknown parameters
u0_guess = [1000.0,750.0,20.0] #initial guess at steady state of state variables, except u0[3] is known to be 20.0
u03 = 20.0 # measured value
plateau = 50.0 # Zero order production rate changes at time = 0, reaches a fractional change given by plateau at start_time
kd = 1.0 # Rate constant for exponential decay back to original value starting at end_time
start_time = 1.0
end_time = 10.0
#### Set up ODE problem
t_end = 36.0 # Length of simulation
tspan = (0.0,t_end) # timespan of simulation
probODE = ODEProblem(sys!,u0_guess,tspan,(p_var,p_fix,kf0,plateau,kd,start_time,end_time))
#### Set up steady state problem
probss = SteadyStateProblem(probODE)
sol_ss = solve(probss)
u0_ss = sol_ss.u
@show "u0_ss true",sol_ss.u
#### Initial guesses
p_guess = [0.1,3.0] #inital guess at values in p_val
kf_guess = [100.0] #initial guess at kf0
termination_condition = NLSolveTerminationCondition(NLSolveTerminationMode.NLSolveDefault;
abstol = 1e-14,
reltol = 1e-14)
#### Set up optimization problem to determine kf for current value of p_guess
set_kf = set_kf_outer(probss,u03,termination_condition)
optkf = OptimizationFunction(set_kf,Optimization.AutoForwardDiff())
prob_set_kf = OptimizationProblem(optkf,kf_guess,(u0_guess,p_guess,p_fix,plateau,kd,start_time,end_time))
@show FiniteDiff.finite_difference_gradient(x -> set_kf(kf_guess,(u0_guess,x,p_fix,plateau,kd,start_time,end_time)),p_guess)
@show ForwardDiff.gradient(x -> set_kf(kf_guess,(u0_guess,x,p_fix,plateau,kd,start_time,end_time)),p_guess )
return nothing
end
main()
Using very small values for abstol and reltol (1e-14) with DynamicSS on a SteadyStateProblem, FiniteDiff solves out to t = 70,000, while ForwardDiff stops at t = 228. The true gradient is [0.0, -296.29629633013286], and FiniteDiff gets that with [-2.6518679901607456e-7, -296.2962963356087], while ForwardDiff is [0.0015469937539040605, -296.29629629628164], good but not great. The gradients are calculated at the bottom of the code with FiniteDiff.finite_difference_gradient and ForwardDiff.gradient
Are the tolerances actually being met with ForwardDiff? Is there a maxiter in the termination conditions in DiffEqBase that is getting hit? Reducing the tolerances to 1e-16 crashes FiniteDiff but ForwardDiff improves to [1.2214211430632144e-5, -296.2962962962962] and goes to t = 1219.
While DynamicSS is not needed on this MWE, it very often is the best choice.