# Setting an upper and lower limit on a differential equation

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?

nope

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!