This is the last piece I need for my model to work, yet the result keeps going over what I thought I had set the limit to.
The simplified scenario is: there is a fishing boat - on the boat work is done - the amount of work depends on the money gained from fishing. There is a limit to how much or how little work is done on the boat. If the upper limit (max_work) is reached a new boat is bought, and if the lower limit (min_work) is reached a boat is sold. (Boat buying/selling is a different equation but I wanted to give more background)
...
if (revenue - cost) > 0.0
limit = max_work
effect = abs(max(0.0, (limit - work)))
else
limit = clearance_min
effect = abs(min(0.0, (limit - work)))
end
du[work] = 365.25*µ*(revenue - cost)*effect
work = u[work]
...
But this doesn’t work, so I tried
if (revenue - cost) > 0.0
limit = max_work
effect = abs(max(0.0, (limit - work)))
else
limit = min_work
effect = abs(min(0.0, (limit - work)))
end
du[work] = 365.25*µ*(revenue - cost)*effect
work = u[work]
if work > max_work
work = max_work
u[work] = max_work
end
if work < min_work
work = min_work
u[work] = min_work
end
And this does not work either. y3 is the upper limit, u84 is work.
See DiffEqCallbacks.jl, specifically the
GeneralDomain(g, u=nothing; nlsolve=NLSOLVEJL_SETUP(), save=true,
abstol=nothing, scalefactor=nothing, autonomous=numargs(g)==2,
nlopts=Dict(:ftol => 10*eps()))
which fixes the solution into a domain. From the docs: “Domains are specified by in-place functions g(u, resid) or g(t, u, resid) that calculate residuals of a state vector u at time t relative to that domain. As for PositiveDomain , steps are accepted if residuals of the extrapolated values at the next time step are below a certain tolerance.”
I think that should do what you are looking for.
Thank you for your reply. What is the reason for not being able to set the limit mathematically?
If your ODE actually goes over the limit, what would this even mean? Do you have an analytical result that says the solution should never go over that limit? The first thing you might want to do is reduce the tolerances a little bit.
1 Like
If I understand your question, the ODE should not be possible to go over the limit. This is due to how the system is designed. I have three base work scenarios: revenue driven, cost minimization, and static, or in-between the other two. A fourth scenario, the one I am trying to plot now, is adaptive based on profit (revenue - cost). It should not be able to go over the system’s maximum work rate.
My solver is defined as
fishing_sol = solve(fishing_prob, alg = CVODE_BDF(), abstol = 1e-6, reltol = 1e-3, callback=cb_set)
I will try reducing the tolerances. Do I need to also state error tolerances in a saving callback?
Reducing tolerances (abstol = 1e-5, reltol = 1e-2) makes it wiggle around near the end, but the issue remains.
No, reduce the numbers: abstol=1e-8,reltol=1e-8.
Right. Sorry, my mistake.
It doesn’t change much I’m afraid.
This is becoming very complicated.
Is this the correct way to define GeneralDomain?
function g(u, resid)
resid[end] = max(0, u[work] - max_work)
end
new_cb = GeneralDomain(g, u=nothing; nlsolve=NLSOLVEJL_SETUP())
And now I need to define NLSOLVEJL_SETUP but it feels this is over my head. https://docs.sciml.ai/release-5.0/features/linear_nonlinear.html#Nonlinear-Solvers:-nlsolve-Specification-1
I still don’t understand why it’s not possible to define the limits mathematically.
Because the ODE has a unique solution given by u' = f(u,p,t). Mathematically, the only way to add limits to the ODE are if the limits are larger than the solution. This means that numerically the limits only make sense to restrict the function to the true limits of the true u(t): you cannot numerically add limits that will change u(t), because then the solver will exit and fail when it realizes that the problem is not well-posed.
This should be correct to about 5-6 decimal places, so if that’s still above where you want, you need to look at the formulation of the model.
Using a callback to change the values gets the required effect…
function condition_cl(bioS, t, integrator)
max_work - bioS[end]
end
function affect_cl!(integrator)
integrator.u[end] = max_work - 0.001
end
…but looks suspicious
It looks like the true solution goes over that “limit”. Are you sure the equation is purely conservative enough to allow for that limit to hold?
1 Like
My colleague figured it out - things were not defined in the right order and stuff was being calculated with the wrong values. 
2 Likes
Good to hear it’s all working now. Good luck!
