ModelingToolkitTolerances.jl is designed to help with the selection of abstol and reltol (see stepsize control for more info) when using ModelingToolkit.jl and DifferentialEquations.jl to solve models. These tolerances act on the error estimates computed when solving the model. These error estimates are very abstract from the physics of the problem and are therefore challenging to set in some cases. As can be seen from the documentation, the defaults for abstol and reltol are 1e-6 and 1e-3, respectively. If these defaults give a bad result, how should they be adjusted? It's very difficult to know what to do in this case, and that is where ModelingToolkitTolerances.jl seeks to help.
Furethermore, in some cases it's sometimes obvious if a model solution was good or bad, but in other cases the question arrises: What is the quality of the model solution? What is it's "error"? With the knowledge that we have from a ModelingToolkit.jl model in addition to the power of Automatic Differentiation, we can now solve for the model error, and even check the error of each solved equation, whether it was algebraic or differential. One of the easiest ways to use this package is to simply check the quality of a model solution. This can be done in a few different ways.
Take the following simple ODE of Euler method with a large step size, we know that this will provide a low quality solution. Let's take a look at the solution with a step size of 2s.
using ModelingToolkit, OrdinaryDiffEq, Plots
using ModelingToolkit: D_nounits as D, t_nounits as t
vars = @variables x(t)=1
eq = D(x) ~ x/4
@mtkcompile sys = System([eq], t, vars, [])
prob = ODEProblem(sys, [], (0, 10))
sol = solve(prob, Euler(); dt=2, adaptive=false)
plot(sol; idxs=x)
Without the analytical solution, we don't know how to quantify the error. With ModelingToolkitTolerances we can now estimate this error. This can be seen visually by adding true (or a number to represent the residual limit) to the 2nd argument of plot(sol::ODESolution, residual_settings::Union{Bool, Real}; kwargs...)
using ModelingToolkitTolerances
plot(sol, true; idxs=x) # residual limit defaults to 1.0
To get the residual information directly we can use the function residual(sol::ODESolution)...
res = residual(sol)
plot(res)
We can convert this to a single number by taking the maximum of the norm.
using LinearAlgebra
r = maximum(norm(res))0.4359462665530558
When using models in an algrithm, say in a loss function for example, it's a good idea to add this to a quality check. For example...
residual_limit = 0.1
accept_result = (sol.retcode == ReturnCode.Success) && (maximum(norm(residual(sol))) < residual_limit)false
If we decrease the step size we should get a better result. Now we see an acceptable result...
sol = solve(prob, Euler(); dt=0.1, adaptive=false)
res = residual(sol)
r = maximum(norm(res)) # 0.05
accept_result = (sol.retcode == ReturnCode.Success) & (r < residual_limit) # true
plot(sol, true; idxs=x)
Take the following hydraulic model of a pressure source, connected to a fixed volume with a valve. The valve starts closed and then opens with a linear ramp. We would expect to see the pressure in the volume rise smoothly, quickly and then slowly, in an S shaped ramp, tabling off at a constant pressure matching the source.
Let's see what we get when we solve this model...
using ModelingToolkit, OrdinaryDiffEq
using ModelingToolkit: t_nounits as t, D_nounits as D
import ModelingToolkitStandardLibrary.Hydraulic.IsothermalCompressible as IC
import ModelingToolkitStandardLibrary.Blocks as B
using Plots
function ValveVolume(; name)
pars = []
systems = @named begin
fluid = IC.HydraulicFluid()
src = IC.FixedPressure(; p = 10e5)
vol = IC.FixedVolume(; vol = 5, p_int=1e5)
valve = IC.Valve(; Cd = 1e5, minimum_area = 0)
ramp = B.Ramp(; height = 0.1, duration = 0.1, offset = 0, start_time = 0.1, smooth = true)
end
eqs = [connect(fluid, src.port)
connect(src.port, valve.port_a)
connect(valve.port_b, vol.port)
connect(valve.area, ramp.output)]
return System(eqs, t, [], pars; name, systems)
end
@mtkcompile sys = ValveVolume()
prob = ODEProblem(sys, [], (0, 5))
sol = solve(prob, Rodas5P()) # <-- using default abstol=1e-6 & reltol=1e-3
plot(sol; idxs=sys.vol.port.p, ylabel="pressure [Pa]")
scatter!(sol.t, sol[sys.vol.port.p]; label="solution points")
What we see here is a very strange solution which seems to be unstable. Clearly we have an invalid solution. Let's use ModelingToolkitTolerances.jl to see what the model error is...
using ModelingToolkitTolerances
resid = residual(sol, 0:0.01:5)
plot(resid)
As can be seen, the error or residual of the model is quite larger than 0. A perfect solution (i.e. an analytical solution) to the problem would give a residual of 0. We can see then that this solution, as expected, is not high quality. So what should we do to get a better solution? Let's use ModelingToolkitTolerances.jl to explore what we can do using the analysis function...
resids = analysis(prob, Rodas5P())Summary of Max Residuals
┌────────┬────────────────┬─────────────────┬─────────────────┬──────────────────┐
│ abstol │ reltol = 0.001 │ reltol = 1.0e-6 │ reltol = 1.0e-9 │ reltol = 1.0e-12 │
├────────┼────────────────┼─────────────────┼─────────────────┼──────────────────┤
│ 0.001 │ 12.2 │ 0.981 │ 0.542 │ 0.543 │
│ 1.0e-6 │ 12.2 │ 1.05 │ 0.0401 │ * 3.52e-5 * │
│ 1.0e-9 │ N/A │ N/A │ N/A │ N/A │
└────────┴────────────────┴─────────────────┴─────────────────┴──────────────────┘
In plotted form...
plot(resids)
What we see here is that we need to adjust the abstol and reltol to 1e-6 and 1e-9 to get the residual to as close to zero as possible. Let's apply that and
see the new result.
sol = solve(prob, Rodas5P(); abstol=1e-6, reltol=1e-9)
Plots.plot(sol; idxs=sys.vol.port.p, ylabel="pressure [Pa]")
Plots.scatter!(sol.t, sol[sys.vol.port.p]; label="solution points")
As can be seen we now get the expected smooth solution to a constant pressure.
We can also get a work precision plot based on the calculated residual.
work_precision(resids)
Here we can see that the point labeled "a6,r12" gives the optimal result, with a minimal residual and a very fast solve time to go with it. Note the label represents the abstol and the reltol exponents, respectively. Therefore, "a6,r12" is translated to "abstol=1e-6, reltol=1e-12".
Often areas of high error will also experience higher compute time. It's also possible to gather this information by adding true as the 2nd positional argument to solve, for example
sol, cpu_timing = solve(prob, true, Rodas5P());
plot(cpu_timing)
Note that we have used @mtkbuild which runs structural_simplify() on the model and reduces the 20 equation system down to 1 equation. So what if we wanted to see the residual of each of those equations individually? This might help us determine which equation is causing the model to struggle. To do this we will use the function no_simplify()...
using ModelingToolkitTolerances: no_simplify
@named sys = ValveVolume()
sys = no_simplify(sys)
st = unknowns(sys)
u0 = st .=> sol(0.0; idxs=st) # <-- create initialization from previous solution
prob = ODEProblem(sys, u0, (0, 1); build_initializeprob=false)
sol = solve(prob, Rodas5P()) # <-- using default abstol & reltol
Plots.plot(sol, 0.05; idxs=sys.vol₊port₊p)
Plots.scatter!(sol.t, sol[sys.vol₊port₊p]; label="solution points")
Here we get a very interesting solution, the adaptive time stepping is working more as expected with default tolerance. Therefore this shows us that structural_simplify() is greatly influencing the solution.
The residual view has also been added to the plot (with the 2nd positional argument for resdiual scale, 0.05). What we see here is the combined residual of all 20 equations. As can be seen the result is much better now, with some increase in residual around the transitions. If we want to look at the algebraic vs. differential equations, we can do the following...
using ModelingToolkitTolerances: ALGEBRAIC, DIFFERENTIAL
resid = residual(sol)
p1 = plot(resid, ALGEBRAIC)
p2 = plot(resid, DIFFERENTIAL)
plot(p1,p2)
As can be seen, all but 1 of the equations are algebraic. And we can see that roughly one equation contributes to the total error, which is equation 9. We can plot specific equations with the following...
using ModelingToolkitTolerances: ResidualSettings
plot(resid, ResidualSettings(9))
We can take a look at the equation...
eqs = equations(sys)
eqs[9]0 ~ -valve₊base₊port_a₊dm(t) + ifelse(valve₊base₊area(t) > valve₊base₊minim
um_area, valve₊base₊area(t), valve₊base₊minimum_area)*ifelse(abs((200.0valv
e₊base₊port_a₊ρ*(valve₊base₊port_a₊p(t) - valve₊base₊port_b₊p(t))) / ifelse
((valve₊base₊port_a₊p(t) - valve₊base₊port_b₊p(t)) > 0, valve₊base₊Cd, valv
e₊base₊Cd_reverse)) >= 1, 0.1(abs((200.0valve₊base₊port_a₊ρ*(valve₊base₊por
t_a₊p(t) - valve₊base₊port_b₊p(t))) / ifelse((valve₊base₊port_a₊p(t) - valv
e₊base₊port_b₊p(t)) > 0, valve₊base₊Cd, valve₊base₊Cd_reverse))^0.5)*sign((
2valve₊base₊port_a₊ρ*(valve₊base₊port_a₊p(t) - valve₊base₊port_b₊p(t))) / i
felse((valve₊base₊port_a₊p(t) - valve₊base₊port_b₊p(t)) > 0, valve₊base₊Cd,
valve₊base₊Cd_reverse)), (20.0valve₊base₊port_a₊ρ*(valve₊base₊port_a₊p(t)
- valve₊base₊port_b₊p(t))) / ifelse((valve₊base₊port_a₊p(t) - valve₊base₊po
rt_b₊p(t)) > 0, valve₊base₊Cd, valve₊base₊Cd_reverse))
This is the valve equation, which is non-linear and makes sense that this is the biggest contributor to the error.
The concept of ModelingToolkitTolerances is quite simple. From ModelingToolkit, we know the symbolic representation of the equations, therefore we know where the derivatives in the system are. Using ForwardDiff and the ODESolution function, we can easily obtain derivatives numerically. Pair those 2 things together and it's easy to calculate the residual of all equations in a system.
Take the simple convection cooling model...
T_inf=300; h=0.7; A=1; m=0.1; c_p=1.2
vars = @variables begin
T(t)=301
end
eqs = [
D(T) ~ (h * A) / (m * c_p) * (T_inf - T)
]
@mtkcompile sys = System(eqs, t, vars, [])
prob = ODEProblem(sys, [], (0, 10))
sol = solve(prob, Tsit5())
plot(sol; idxs=T)
We are plotting the temperature, which should be smoothly cooling from 301K to 300K and reaching steady state. But something in the numerical solution goes wrong and the temperature becomes unstable. In this case we clearly have a measurable residual error in the solved equation. The residual is calculated by simply taking the left hand side of the equation subtracted from the right hand side:
res = ( (h * A)/(m * c_p) * (T_inf - T) ) - ( D(T) )From the numerical solution, we know T, but we don't have D(T). As mentioned, we can estimate this using ForwardDiff...
using ForwardDiff
Tsol(t) = sol(t; idxs=sys.T)
dTsol(t) = ForwardDiff.derivative(Tsol, t)We can now calculate the residual as follows...
resf(t) = ( (h * A) / (m * c_p) * (T_inf - Tsol(t)) ) - ( dTsol(t) )
times = 0:1e-3:10
plot(times, resf; label="manual")
This is exactly what is computed automatically when asking for the residual from ModelingToolkitTolerances
res = residual(sol, times)
plot!(res.t, res.residuals[:,1]; label="ModelingToolkitTolerances")
Note: calling plot(info::ResidualInfo, settings::ResidualSettings = SUMMARY) will plot the norm of all requested equations in settings.
Using the analysis function we can quickly find a better tolerance setting...
resids = analysis(prob, Tsit5())Summary of Max Residuals
┌────────┬────────────────┬─────────────────┬─────────────────┬──────────────────┐
│ abstol │ reltol = 0.001 │ reltol = 1.0e-6 │ reltol = 1.0e-9 │ reltol = 1.0e-12 │
├────────┼────────────────┼─────────────────┼─────────────────┼──────────────────┤
│ 0.001 │ 2.6 │ 0.0104 │ 0.0068 │ 0.00678 │
│ 1.0e-6 │ 2.58 │ 0.00224 │ 9.14e-6 │ 7.32e-6 │
│ 1.0e-9 │ 2.58 │ 0.00227 │ 2.17e-6 │ * 2.01e-8 * │
└────────┴────────────────┴─────────────────┴─────────────────┴──────────────────┘
plot(resids)
As can be seen, to get a residual less than 1 simply requires setting the reltol to 1e-6.
sol = solve(prob, Tsit5(); reltol=1e-6)
plot(sol; idxs=T)
As can be seen, we now have a stable solution, as expected.
Using saveat can affect the quality of the ForwardDiff.derivative computed from the ODESolution. Avoid setting this keyword when using residual() or analysis().
residual(sol::ODESolution, tms = default_range(sol); abstol=0.0, reltol=0.0, timing=0.0)
Calculates residual information for sol::ODESolution coming from a ModelingToolkit model at the times tms, which defaults to the solved time points sol.t. Returns as ResidualInfo object.
Keyword Arguments (used by analysis() function):
abstol: this simply records theabstolthat was used to createsolobject. Since this information is not recorded in anODESolutionobject it must be specified explicitly for reference. Used byanalysis()function.reltol: same as above, but forreltoltiming: this records the corresponding solution time reference. Used byanalysis()function.
analysis(prob::ODEProblem, solver, tms = collect(prob.tspan[1]:1e-3:prob.tspan[2]); kwargs...)
Runs a 3 x 3 study of abstol and reltol = [1e-3, 1e-6, 1e-9]. Returns a Vector{ResidualInfo} which can be sent to plot() and work_precision() functions for visual analysis. A solver from OrdinaryDiffEq must be passed as the 2nd argument. The time points used for residual calculation are provided by tms and default to the prob.tspan spaced by 1ms. Provided kwargs... are passed to the solve