# 7. Rigorous error bounds using IntervalArithmetic

This example is also available as a Jupyter notebook that can be run locally The notebook can be found in the `examples`

directory of the package. If the notebooks are missing, you may need to run `using Pkg; Pkg.build()`

.

## Floating-point error

In computers, real numbers are commonly approximated using floating-point numbers, such as Julia's `Float64`

. Unfortunately, not all real numbers can be exactly represented as a finite-size floating-point number, and the results of operations on floating-point numbers can only approximate the results of applying the operation to a true real number. This results in peculiarities like:

`2.6 - 0.7 - 1.9`

2.220446049250313e-16

IntervalArithmetic.jl can be used to quantify floating point error, by computing *rigorous* worst-case bounds on floating point error, within which the true result is *guaranteed* to lie.

`using IntervalArithmetic`

Activating environment at `~/work/RigidBodyDynamics.jl/RigidBodyDynamics.jl/examples/7. Rigorous error bounds using IntervalArithmetic/Project.toml`

IntervalArithmetic.jl provides the `Interval`

type, which stores an upper and a lower bound:

`i = Interval(1.0, 2.0)`

[1, 2]

`dump(i)`

IntervalArithmetic.Interval{Float64} lo: Float64 1.0 hi: Float64 2.0

IntervalArithmetic.jl provides overloads for most common Julia functions that take these bounds into account. For example:

`i + i`

[2, 4]

`sin(i)`

[0.84147, 1]

Note that the bounds computed by IntervalArithmetic.jl take floating point error into account. Also note that a given real number, once converted to (approximated by) a floating-point number may not be equal to the original real number. To rigorously construct an `Interval`

that contains a given real number as an input, IntervalArithmetic.jl provides the `@interval`

macro:

```
i = @interval(2.9)
i.lo === i.hi
```

false

`dump(i)`

IntervalArithmetic.Interval{Float64} lo: Float64 2.9 hi: Float64 2.9000000000000004

Compare this to

```
i = Interval(2.9)
i.lo === i.hi
```

true

`dump(i)`

IntervalArithmetic.Interval{Float64} lo: Float64 2.9 hi: Float64 2.9

As an example, consider again the peculiar result from before, now using interval arithmetic:

`i = @interval(2.6) - @interval(0.7) - @interval(1.9)`

[-6.66134e-16, 2.22045e-16]

showing that the true result, `0`

, is indeed in the guaranteed interval, and indeed:

```
using Test
@test (2.6 - 0.7 - 1.9) ∈ i
```

Test Passed

## Accuracy of RigidBodyDynamics.jl's `mass_matrix`

Let's use IntervalArithmetic.jl to establish rigorous bounds on the accuracy of the accuracy of the `mass_matrix`

algorithm for the Acrobot (double pendulum) in a certain configuration. Let's get started.

`using RigidBodyDynamics`

We'll create a `Mechanism`

by parsing the Acrobot URDF, passing in `Interval{Float64}`

as the type used to store the parameters (inertias, link lengths, etc.) of the mechanism. Note that the parameters parsed from the URDF are treated as floating point numbers (i.e., like `Interval(2.9)`

instead of `@interval(2.9)`

above).

```
const T = Interval{Float64}
srcdir = dirname(pathof(RigidBodyDynamics))
urdf = joinpath(srcdir, "..", "test", "urdf", "Acrobot.urdf")
const mechanism = parse_urdf(urdf; scalar_type=T)
state = MechanismState(mechanism)
```

MechanismState{IntervalArithmetic.Interval{Float64}, IntervalArithmetic.Interval{Float64}, IntervalArithmetic.Interval{Float64}, …}(…)

Let's set the initial joint angle of the shoulder joint to the smallest `Interval{Float64}`

containing the real number $1$, and similarly for the elbow joint:

```
shoulder, elbow = joints(mechanism)
set_configuration!(state, shoulder, @interval(1))
set_configuration!(state, elbow, @interval(2));
```

And now we can compute the mass matrix as normal:

`M = mass_matrix(state)`

2×2 LinearAlgebra.Symmetric{IntervalArithmetic.Interval{Float64},Array{IntervalArithmetic.Interval{Float64},2}}: [1.8307, 1.83071] [0.913853, 0.913854] [0.913853, 0.913854] [1.32999, 1.33001]

Woah, those bounds look pretty big. RigidBodyDynamics.jl must not be very accurate! Actually, things aren't so bad; the issue is just that IntervalArithmetic.jl isn't kidding when it comes to guaranteed bounds, and that includes printing the numbers in shortened form. Here are the lengths of the intervals:

`err = map(x -> x.hi - x.lo, M)`

2×2 Array{Float64,2}: 3.9968e-15 4.88498e-15 4.88498e-15 6.43929e-15

`@test maximum(abs, err) ≈ 0 atol = 1e-14`

Test Passed

## Rigorous (worst-case) uncertainty propagation

IntervalArithmetic.jl can also be applied to propagate uncertainty in a rigorous way when the inputs themselves are uncertain. Consider for example the case that we only know the joint angles up to $\pm 0.05$ radians:

```
set_configuration!(state, shoulder, @interval(0.95, 1.05))
set_configuration!(state, elbow, @interval(1.95, 2.05));
```

and let's compute bounds on the center of mass position:

`center_of_mass(state)`

Point3D in "world": IntervalArithmetic.Interval{Float64}[[-0.770193, -0.630851], [0.199999, 0.200001], [0.0167304, 0.163822]]

Note that the bounds on the $y$-coordinate are very tight, since our mechanism only lives in the $x$-$z$ plane.

*This page was generated using Literate.jl.*